{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2f6cd6d0-1463-43ab-aca0-72f6dd184017",
   "metadata": {},
   "outputs": [],
   "source": [
    "# options(repos = c(CRAN = \"https://mirrors.aliyun.com/CRAN/\"))\n",
    "# install.packages(\"coin\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "76195f31-2574-4e67-8dc7-2aa55e92eab0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   treatment score\n",
      "1          A    40\n",
      "2          A    57\n",
      "3          A    45\n",
      "4          A    55\n",
      "5          A    58\n",
      "6          B    57\n",
      "7          B    64\n",
      "8          B    55\n",
      "9          B    62\n",
      "10         B    65\n"
     ]
    }
   ],
   "source": [
    "score <- c(40,57,45,55,58,57,64,55,62,65)\n",
    "treatment <- factor(c(rep(\"A\",5),rep(\"B\",5)))\n",
    "mydata <- data.frame(treatment,score)\n",
    "print(mydata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "434ae4c0-b938-414b-8226-fee7ef39437a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "\tTwo Sample t-test\n",
       "\n",
       "data:  score by treatment\n",
       "t = -2.345, df = 8, p-value = 0.04705\n",
       "alternative hypothesis: true difference in means between group A and group B is not equal to 0\n",
       "95 percent confidence interval:\n",
       " -19.0405455  -0.1594545\n",
       "sample estimates:\n",
       "mean in group A mean in group B \n",
       "           51.0            60.6 \n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "t.test(score~treatment,data=mydata,var.equal=TRUE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "10744df1-8034-4a59-8616-5d6ec338eb84",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "\"package 'coin' was built under R version 4.4.3\"\n",
      "Loading required package: survival\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "\tExact Two-Sample Fisher-Pitman Permutation Test\n",
       "\n",
       "data:  score by treatment (A, B)\n",
       "Z = -1.9147, p-value = 0.07143\n",
       "alternative hypothesis: true mu is not equal to 0\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "library(coin)\n",
    "library(survival)\n",
    "oneway_test(score~treatment,data=mydata,distribution=\"exact\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "1d6af3ac-234e-417b-bd1f-8e700fbe0595",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "\tWilcoxon rank sum exact test\n",
       "\n",
       "data:  Prob by So\n",
       "W = 81, p-value = 8.488e-05\n",
       "alternative hypothesis: true location shift is not equal to 0\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "library(MASS)\n",
    "\n",
    "UScrime <- transform(UScrime,So=factor(So))\n",
    "wilcox.test(Prob~So,data=UScrime,distribution=\"exact\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ec15c1d8-1a31-48a7-9dd0-04967e96e4f0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "\tExact Wilcoxon-Mann-Whitney Test\n",
       "\n",
       "data:  Prob by So (0, 1)\n",
       "Z = -3.7493, p-value = 8.488e-05\n",
       "alternative hypothesis: true mu is not equal to 0\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "wilcox_test(Prob~So,data=UScrime,distribution=\"exact\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "87a4b6d0-63b0-456a-b3ff-2ff63be6fd24",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "16"
      ],
      "text/latex": [
       "16"
      ],
      "text/markdown": [
       "16"
      ],
      "text/plain": [
       "[1] 16"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nrow(UScrime[UScrime$So==1,])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a6498f5-e529-47ad-8ab2-529efc65d44e",
   "metadata": {},
   "source": [
    "按照公式计算的结果和显示的结果不一样？"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "e8194559-fe09-436a-bdd6-64a9b7e22388",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "-6.8025508142641"
      ],
      "text/latex": [
       "-6.8025508142641"
      ],
      "text/markdown": [
       "-6.8025508142641"
      ],
      "text/plain": [
       "[1] -6.802551"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "(81 - 16*48/2)/sqrt(16*31*48/12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "99ec4c67-ef79-42aa-97f4-14d025da5373",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Loading required package: mvtnorm\n",
      "\n",
      "Loading required package: TH.data\n",
      "\n",
      "\n",
      "Attaching package: 'TH.data'\n",
      "\n",
      "\n",
      "The following object is masked from 'package:MASS':\n",
      "\n",
      "    geyser\n",
      "\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "\tApproximative K-Sample Fisher-Pitman Permutation Test\n",
       "\n",
       "data:  response by\n",
       "\t trt (1time, 2times, 4times, drugD, drugE)\n",
       "chi-squared = 36.381, p-value < 1e-04\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "library(multcomp)\n",
    "set.seed(1234)\n",
    "oneway_test(response~trt,data=cholesterol,distribution=approximate(nresample=9999))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f42c4f42-9022-4556-a6d1-be5b91058301",
   "metadata": {},
   "source": [
    "到目前为止，给我感觉，置换检验更像是传统检验的一种抽样估计……"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "992eaf3c-c27e-42c7-8f17-14c98d25baf4",
   "metadata": {},
   "source": [
    "\n",
    "置换检验应用于卡方检验（通常用于列联表数据，检验两个分类变量的独立性）时，核心思想是通过随机重排数据来模拟原假设成立时的分布，从而避免依赖传统卡方检验的理论分布假设（如大样本下的卡方近似）。以下是具体流程：\n",
    "\n",
    "\n",
    "**1. 明确原假设与备择假设**\n",
    "\n",
    "- **原假设 $H_0$**：两个分类变量相互独立（即无关联）。  \n",
    "- **备择假设 $H_1$**：两个分类变量存在关联（单侧或双侧检验，根据需求确定）。\n",
    "\n",
    "\n",
    "**2. 构建列联表并计算观测卡方统计量**\n",
    "\n",
    "- **观测列联表**：设数据为一个 $R \\times C$ 的列联表，记录两个分类变量的频数，记为 $O_{ij}$（第 $i$ 行、第 $j$ 列的观测频数）。  \n",
    "- **计算期望频数**：在原假设（独立）下，期望频数 $E_{ij} = \\frac{\\text{第 }i\\text{ 行合计} \\times \\text{第 }j\\text{ 列合计}}{\\text{总样本量 }n}$。  \n",
    "- **观测卡方统计量**：  \n",
    "  $$\n",
    "  \\chi^2_{\\text{obs}} = \\sum_{i=1}^R \\sum_{j=1}^C \\frac{(O_{ij} - E_{ij})^2}{E_{ij}}\n",
    " $$\n",
    "\n",
    "\n",
    "**3. 置换重排数据以模拟原假设**\n",
    "\n",
    "在原假设（变量独立）下，数据的行或列标签可随机置换，以保持边际分布（行和与列和）不变，模拟变量无关联的情况。具体步骤如下：\n",
    "\n",
    "**置换方法（以两分类变量为例）**  \n",
    "\n",
    "- **场景**：假设其中一个变量是分组变量（如处理组/对照组），另一个是结果变量（如成功/失败），形成 $2 \\times 2$ 列联表。  \n",
    "- **置换操作**：  \n",
    "  1. 固定所有样本的结果变量（列标签），随机重排分组变量（行标签），即将每个样本的分组标签（如处理/对照）随机分配，不改变结果变量的分布。  \n",
    "  2. 对于更一般的 $R \\times C$ 表，可随机置换其中一个变量的分类标签（如行标签），确保行和与列和不变（即保持边际频数固定），从而维持原假设下的独立性。  \n",
    "\n",
    "**重复置换** \n",
    "\n",
    "- 进行 $B$ 次独立置换（如 $B = 10^4$ 次），每次置换后生成新的列联表，记为 $O^{(b)}_{ij}$（$b = 1, 2, \\dots, B$）。  \n",
    "- 对每次置换后的列联表，计算卡方统计量 $\\chi^2_{(b)}$。\n",
    "\n",
    "\n",
    "**4. 构建置换分布并计算 $p$ 值**  \n",
    "\n",
    "- **置换分布**：收集 $B$ 次置换得到的卡方统计量 $\\{\\chi^2_{(1)}, \\chi^2_{(2)}, \\dots, \\chi^2_{(B)}\\}$，形成经验分布。  \n",
    "- **计算 $p$ 值**：  \n",
    "  - 若为双侧检验（关联方向不确定）：  \n",
    "    $$\n",
    "    p = \\frac{\\text{置换中满足 } \\chi^2_{(b)} \\geq \\chi^2_{\\text{obs}} \\text{ 的次数} + 1}{B + 1}\n",
    "    $$ \n",
    "    （加1是为了避免 $p=0$，考虑观测值本身是否计入）  \n",
    "  - 若为单侧检验（如正向关联）：根据备择假设方向调整（仅计算大于或小于观测值的情况）。  \n",
    "\n",
    "\n",
    "**5. 决策与结论**  \n",
    "\n",
    "- 比较 $p$ 值与显著性水平 $\\alpha$（如 0.05）：  \n",
    "  - 若 $p \\leq \\alpha$：拒绝原假设，认为两变量存在显著关联。  \n",
    "  - 若 $p > \\alpha$：不拒绝原假设，无足够证据表明两变量关联。  \n",
    "\n",
    "\n",
    "\n",
    "**关键特点与注意事项**  \n",
    "\n",
    "1. **优势**：  \n",
    "   - 不依赖卡方分布的大样本近似，适用于小样本或期望频数较小的场景（传统卡方检验要求期望频数 ≥ 5，否则结果不可靠）。  \n",
    "   - 精确计算 $p$ 值，基于数据本身的置换分布，而非理论分布。  \n",
    "\n",
    "2. **置换逻辑**：  \n",
    "   - 置换需保持边际频数不变（行和与列和固定），否则会破坏原假设下的独立性前提。  \n",
    "   - 计算量随 $B$ 增大而增加，但现代计算能力可处理较大 $B$（如 $10^5$ 次）。  \n",
    "\n",
    "3. **与传统卡方检验的区别**：  \n",
    "   - 传统方法使用理论卡方分布近似 $p$ 值，置换检验使用经验置换分布精确计算 $p$ 值（尤其在小样本时更可靠）。  \n",
    "\n",
    "\n",
    "**示例：两独立样本的 $2 \\times 2$ 表置换检验**  \n",
    "\n",
    "- 观测数据：处理组和对照组的成功/失败频数。  \n",
    "- 置换操作：固定成功/失败总数，随机重排每个样本属于处理组或对照组，生成新的列联表，计算卡方统计量。  \n",
    "- 通过 $B$ 次置换，统计有多少次置换卡方 ≥ 观测卡方，得到 $p$ 值。  \n",
    "\n",
    "通过以上流程，置换检验为卡方检验提供了一种不依赖分布假设的精确推断方法，尤其适用于传统方法不满足条件的场景。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "22aa62cd-403e-42c0-af81-58654713f79f",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Loading required package: grid\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "\tApproximative Pearson Chi-Squared Test\n",
       "\n",
       "data:  Treatment by Improved (1, 2, 3)\n",
       "chi-squared = 13.055, p-value = 0.0018\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "library(coin)\n",
    "library(vcd)\n",
    "set.seed(1234)\n",
    "Arthritis <- transform(Arthritis,Improved=as.factor(as.numeric(Improved)))\n",
    "chisq_test(Treatment~Improved,data=Arthritis,distribution=approximate(nresample=9999))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "8fc6e32d-15bb-4297-b98c-ae34ea7d2e1b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "\tApproximative Spearman Correlation Test\n",
       "\n",
       "data:  Illiteracy by Murder\n",
       "Z = 4.7065, p-value < 1e-04\n",
       "alternative hypothesis: true rho is not equal to 0\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "states <- as.data.frame(state.x77)\n",
    "set.seed(1234)\n",
    "coin::spearman_test(Illiteracy~Murder,data=states,distribution=approximate(nresample=9999))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "c661e729-ca2b-4b57-b0b3-30873f477b51",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "1.26002992772989e-06"
      ],
      "text/latex": [
       "1.26002992772989e-06"
      ],
      "text/markdown": [
       "1.26002992772989e-06"
      ],
      "text/plain": [
       "[1] 1.26003e-06"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "1-pnorm(abs(4.7065))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "2b7249a9-3194-4f22-a178-f7dcdaff5f1a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    U1 U2\n",
      "1  108 41\n",
      "2   96 36\n",
      "3   94 33\n",
      "4  102 39\n",
      "5   91 20\n",
      "6   84 29\n",
      "7   97 38\n",
      "8   79 35\n",
      "9   81 28\n",
      "10 100 24\n"
     ]
    }
   ],
   "source": [
    "library(MASS)\n",
    "\n",
    "print(UScrime[1:10,c(\"U1\",\"U2\")])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "aa77e551-bfe8-4739-97ca-2526909ca14f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "\tExact Wilcoxon-Pratt Signed-Rank Test\n",
       "\n",
       "data:  y by x (pos, neg) \n",
       "\t stratified by block\n",
       "Z = 5.9691, p-value = 1.421e-14\n",
       "alternative hypothesis: true mu is not equal to 0\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "wilcoxsign_test(U1~U2,data=UScrime,distribution=\"exact\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d2aad6d-0e59-482d-8e3d-566816f0502f",
   "metadata": {},
   "source": [
    "Wilcoxon 符号秩检验主要用于配对样本的非参数检验，用于判断配对样本之间是否存在显著差异。其秩的计算步骤如下：\n",
    "\n",
    "**步骤 1：计算配对样本的差值**\n",
    "\n",
    "假设有配对样本$(X_1,Y_1),(X_2,Y_2),\\cdots,(X_n,Y_n)$，首先计算每对样本的差值$d_i = X_i - Y_i$，$i = 1,2,\\cdots,n$。\n",
    "\n",
    "**步骤 2：去除差值为 0 的样本对**\n",
    "\n",
    "将差值$d_i = 0$的样本对去除，剩下的样本对数量记为$n'$。因为差值为 0 对于判断两组数据是否有差异没有提供有效信息。\n",
    "\n",
    "**步骤 3：计算差值的绝对值**\n",
    "\n",
    "对剩下的$n'$个差值$d_i$，计算它们的绝对值$\\vert d_i\\vert$。\n",
    "\n",
    "**步骤 4：对绝对值进行排序并赋予秩次**\n",
    "\n",
    "将$\\vert d_i\\vert$从小到大进行排序，然后为每个$\\vert d_i\\vert$赋予一个秩次。秩次就是该绝对值在排序后的序列中所处的位置。\n",
    "| 样本对编号 | $X$| $Y$| $d = X - Y$| $\\vert d\\vert$| 排序后的$\\vert d\\vert$| 秩次 | 带符号的秩次 |\n",
    "| ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- |\n",
    "| 1 | 12 | 10 | 2 | 2 | 4 | 4.5 | 4.5 |\n",
    "| 2 | 15 | 16 | -1 | 1 | 1 | 2 | -2 |\n",
    "| 3 | 18 | 17 | 1 | 1 | 1 | 2 | 2 |\n",
    "| 4 | 20 | 22 | -2 | 2 | 4 | 4.5 | -4.5 |\n",
    "| 5 | 22 | 21 | 1 | 1 | 1 | 2 | 2 |\n",
    "\n",
    "**步骤 5：恢复差值的符号**\n",
    "\n",
    "将步骤 4 中得到的秩次乘以对应的差值$d_i$的符号（正差值对应的秩次为正，负差值对应的秩次为负），得到带符号的秩次。\n",
    "\n",
    "**示例计算**\n",
    "\n",
    "假设有以下配对样本数据：\n",
    "| 样本对编号 | $X$| $Y$| $d = X - Y$| $\\vert d\\vert$| 排序后的 $\\vert d\\vert$| 秩次 | 带符号的秩次 |\n",
    "| ---- | ---- | ---- | ----: | ----:| ----: | ----: | ----: |\n",
    "| 1 | 12 | 10 | 2 | 2 | 1 | 1 | 1 |\n",
    "| 2 | 15 | 16 | -1 | 1 | 2 | 2 | -2 |\n",
    "| 3 | 18 | 17 | 1 | 1 | 2 | 2 | 2 |\n",
    "| 4 | 20 | 22 | -2 | 2 | 4 | 4 | -4 |\n",
    "| 5 | 22 | 21 | 1 | 1 | 2 | 2 | 2 |\n",
    "\n",
    "在这个例子中，有三个$\\vert d\\vert$的值为 1，它们原本的秩次应该是 1、2、3，所以它们的共同秩次为$(1 + 2 + 3)\\div3 = 2$；有两个$\\vert d\\vert$的值为 2，它们原本的秩次应该是 4、5，所以它们的共同秩次为$(4 + 5)\\div2 = 4.5$。\n",
    "\n",
    "最后，根据差值的符号得到带符号的秩次。在 Wilcoxon 符号秩检验中，会基于这些带符号的秩次进一步计算检验统计量，以判断配对样本之间是否存在显著差异。 "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "caf4ff11-59ab-4eac-843d-465cca291fe5",
   "metadata": {},
   "source": [
    "### **lmPerm**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41ef7924-fd8a-4c76-8940-b41bd55a26a4",
   "metadata": {},
   "source": [
    "$$F=\\frac{MSR}{MSE} = \\frac{\\frac{SSR}{k}}{\\frac{SSE}{n-k-1}}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e473653d-226e-4a1e-b29d-6c53cec96efc",
   "metadata": {},
   "source": [
    "$SSR = \\sum\\limits_{i=1}^{n}(\\hat y_i - \\bar y)^2$\n",
    "\n",
    "$SSE = \\sum\\limits_{i=1}^{n}( y_i - \\hat y)^2$\n",
    "\n",
    "ps：这里要学习的权重为:$\\beta_0,\\beta_1,\\cdots,\\beta_k$，共有$k+1$个。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "6bcc4d6e-69ad-4317-83fe-8f5ff29bbbdb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "Call:\n",
       "lm(formula = weight ~ height, data = women)\n",
       "\n",
       "Residuals:\n",
       "    Min      1Q  Median      3Q     Max \n",
       "-1.7333 -1.1333 -0.3833  0.7417  3.1167 \n",
       "\n",
       "Coefficients:\n",
       "             Estimate Std. Error t value Pr(>|t|)    \n",
       "(Intercept) -87.51667    5.93694  -14.74 1.71e-09 ***\n",
       "height        3.45000    0.09114   37.85 1.09e-14 ***\n",
       "---\n",
       "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n",
       "\n",
       "Residual standard error: 1.525 on 13 degrees of freedom\n",
       "Multiple R-squared:  0.991,\tAdjusted R-squared:  0.9903 \n",
       "F-statistic:  1433 on 1 and 13 DF,  p-value: 1.091e-14\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fit = lm(weight~height,data=women)\n",
    "summary(fit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "7365165b-eb33-4b5b-802c-9c5d6f4637ae",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "Call:\n",
       "lm(formula = weight ~ height + I(height^2), data = women)\n",
       "\n",
       "Residuals:\n",
       "     Min       1Q   Median       3Q      Max \n",
       "-0.50941 -0.29611 -0.00941  0.28615  0.59706 \n",
       "\n",
       "Coefficients:\n",
       "             Estimate Std. Error t value Pr(>|t|)    \n",
       "(Intercept) 261.87818   25.19677  10.393 2.36e-07 ***\n",
       "height       -7.34832    0.77769  -9.449 6.58e-07 ***\n",
       "I(height^2)   0.08306    0.00598  13.891 9.32e-09 ***\n",
       "---\n",
       "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n",
       "\n",
       "Residual standard error: 0.3841 on 12 degrees of freedom\n",
       "Multiple R-squared:  0.9995,\tAdjusted R-squared:  0.9994 \n",
       "F-statistic: 1.139e+04 on 2 and 12 DF,  p-value: < 2.2e-16\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fit = lm(weight~height+I(height^2),data=women)\n",
    "summary(fit)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f150b05-0c75-4782-af20-44ff087cef15",
   "metadata": {},
   "source": [
    "**所以安装lmPerm的意义在哪里？结果完全一样！！！**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "fffddb91-0bc3-4223-8836-568bcb131695",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"dataframe\">\n",
       "<caption>A data.frame: 10 × 8</caption>\n",
       "<thead>\n",
       "\t<tr><th></th><th scope=col>Population</th><th scope=col>Income</th><th scope=col>Illiteracy</th><th scope=col>Life Exp</th><th scope=col>Murder</th><th scope=col>HS Grad</th><th scope=col>Frost</th><th scope=col>Area</th></tr>\n",
       "\t<tr><th></th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>Alabama</th><td> 3615</td><td>3624</td><td>2.1</td><td>69.05</td><td>15.1</td><td>41.3</td><td> 20</td><td> 50708</td></tr>\n",
       "\t<tr><th scope=row>Alaska</th><td>  365</td><td>6315</td><td>1.5</td><td>69.31</td><td>11.3</td><td>66.7</td><td>152</td><td>566432</td></tr>\n",
       "\t<tr><th scope=row>Arizona</th><td> 2212</td><td>4530</td><td>1.8</td><td>70.55</td><td> 7.8</td><td>58.1</td><td> 15</td><td>113417</td></tr>\n",
       "\t<tr><th scope=row>Arkansas</th><td> 2110</td><td>3378</td><td>1.9</td><td>70.66</td><td>10.1</td><td>39.9</td><td> 65</td><td> 51945</td></tr>\n",
       "\t<tr><th scope=row>California</th><td>21198</td><td>5114</td><td>1.1</td><td>71.71</td><td>10.3</td><td>62.6</td><td> 20</td><td>156361</td></tr>\n",
       "\t<tr><th scope=row>Colorado</th><td> 2541</td><td>4884</td><td>0.7</td><td>72.06</td><td> 6.8</td><td>63.9</td><td>166</td><td>103766</td></tr>\n",
       "\t<tr><th scope=row>Connecticut</th><td> 3100</td><td>5348</td><td>1.1</td><td>72.48</td><td> 3.1</td><td>56.0</td><td>139</td><td>  4862</td></tr>\n",
       "\t<tr><th scope=row>Delaware</th><td>  579</td><td>4809</td><td>0.9</td><td>70.06</td><td> 6.2</td><td>54.6</td><td>103</td><td>  1982</td></tr>\n",
       "\t<tr><th scope=row>Florida</th><td> 8277</td><td>4815</td><td>1.3</td><td>70.66</td><td>10.7</td><td>52.6</td><td> 11</td><td> 54090</td></tr>\n",
       "\t<tr><th scope=row>Georgia</th><td> 4931</td><td>4091</td><td>2.0</td><td>68.54</td><td>13.9</td><td>40.6</td><td> 60</td><td> 58073</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A data.frame: 10 × 8\n",
       "\\begin{tabular}{r|llllllll}\n",
       "  & Population & Income & Illiteracy & Life Exp & Murder & HS Grad & Frost & Area\\\\\n",
       "  & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl>\\\\\n",
       "\\hline\n",
       "\tAlabama &  3615 & 3624 & 2.1 & 69.05 & 15.1 & 41.3 &  20 &  50708\\\\\n",
       "\tAlaska &   365 & 6315 & 1.5 & 69.31 & 11.3 & 66.7 & 152 & 566432\\\\\n",
       "\tArizona &  2212 & 4530 & 1.8 & 70.55 &  7.8 & 58.1 &  15 & 113417\\\\\n",
       "\tArkansas &  2110 & 3378 & 1.9 & 70.66 & 10.1 & 39.9 &  65 &  51945\\\\\n",
       "\tCalifornia & 21198 & 5114 & 1.1 & 71.71 & 10.3 & 62.6 &  20 & 156361\\\\\n",
       "\tColorado &  2541 & 4884 & 0.7 & 72.06 &  6.8 & 63.9 & 166 & 103766\\\\\n",
       "\tConnecticut &  3100 & 5348 & 1.1 & 72.48 &  3.1 & 56.0 & 139 &   4862\\\\\n",
       "\tDelaware &   579 & 4809 & 0.9 & 70.06 &  6.2 & 54.6 & 103 &   1982\\\\\n",
       "\tFlorida &  8277 & 4815 & 1.3 & 70.66 & 10.7 & 52.6 &  11 &  54090\\\\\n",
       "\tGeorgia &  4931 & 4091 & 2.0 & 68.54 & 13.9 & 40.6 &  60 &  58073\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A data.frame: 10 × 8\n",
       "\n",
       "| <!--/--> | Population &lt;dbl&gt; | Income &lt;dbl&gt; | Illiteracy &lt;dbl&gt; | Life Exp &lt;dbl&gt; | Murder &lt;dbl&gt; | HS Grad &lt;dbl&gt; | Frost &lt;dbl&gt; | Area &lt;dbl&gt; |\n",
       "|---|---|---|---|---|---|---|---|---|\n",
       "| Alabama |  3615 | 3624 | 2.1 | 69.05 | 15.1 | 41.3 |  20 |  50708 |\n",
       "| Alaska |   365 | 6315 | 1.5 | 69.31 | 11.3 | 66.7 | 152 | 566432 |\n",
       "| Arizona |  2212 | 4530 | 1.8 | 70.55 |  7.8 | 58.1 |  15 | 113417 |\n",
       "| Arkansas |  2110 | 3378 | 1.9 | 70.66 | 10.1 | 39.9 |  65 |  51945 |\n",
       "| California | 21198 | 5114 | 1.1 | 71.71 | 10.3 | 62.6 |  20 | 156361 |\n",
       "| Colorado |  2541 | 4884 | 0.7 | 72.06 |  6.8 | 63.9 | 166 | 103766 |\n",
       "| Connecticut |  3100 | 5348 | 1.1 | 72.48 |  3.1 | 56.0 | 139 |   4862 |\n",
       "| Delaware |   579 | 4809 | 0.9 | 70.06 |  6.2 | 54.6 | 103 |   1982 |\n",
       "| Florida |  8277 | 4815 | 1.3 | 70.66 | 10.7 | 52.6 |  11 |  54090 |\n",
       "| Georgia |  4931 | 4091 | 2.0 | 68.54 | 13.9 | 40.6 |  60 |  58073 |\n",
       "\n"
      ],
      "text/plain": [
       "            Population Income Illiteracy Life Exp Murder HS Grad Frost Area  \n",
       "Alabama      3615      3624   2.1        69.05    15.1   41.3     20    50708\n",
       "Alaska        365      6315   1.5        69.31    11.3   66.7    152   566432\n",
       "Arizona      2212      4530   1.8        70.55     7.8   58.1     15   113417\n",
       "Arkansas     2110      3378   1.9        70.66    10.1   39.9     65    51945\n",
       "California  21198      5114   1.1        71.71    10.3   62.6     20   156361\n",
       "Colorado     2541      4884   0.7        72.06     6.8   63.9    166   103766\n",
       "Connecticut  3100      5348   1.1        72.48     3.1   56.0    139     4862\n",
       "Delaware      579      4809   0.9        70.06     6.2   54.6    103     1982\n",
       "Florida      8277      4815   1.3        70.66    10.7   52.6     11    54090\n",
       "Georgia      4931      4091   2.0        68.54    13.9   40.6     60    58073"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "head(states,10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "58f01eed-6205-4554-b10b-6f7fd20617f8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "Call:\n",
       "lm(formula = Murder ~ Population + Illiteracy + Income + Frost, \n",
       "    data = states)\n",
       "\n",
       "Residuals:\n",
       "    Min      1Q  Median      3Q     Max \n",
       "-4.7960 -1.6495 -0.0811  1.4815  7.6210 \n",
       "\n",
       "Coefficients:\n",
       "             Estimate Std. Error t value Pr(>|t|)    \n",
       "(Intercept) 1.235e+00  3.866e+00   0.319   0.7510    \n",
       "Population  2.237e-04  9.052e-05   2.471   0.0173 *  \n",
       "Illiteracy  4.143e+00  8.744e-01   4.738 2.19e-05 ***\n",
       "Income      6.442e-05  6.837e-04   0.094   0.9253    \n",
       "Frost       5.813e-04  1.005e-02   0.058   0.9541    \n",
       "---\n",
       "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n",
       "\n",
       "Residual standard error: 2.535 on 45 degrees of freedom\n",
       "Multiple R-squared:  0.567,\tAdjusted R-squared:  0.5285 \n",
       "F-statistic: 14.73 on 4 and 45 DF,  p-value: 9.133e-08\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fit = lm(Murder~Population+Illiteracy+Income+Frost,data=states)\n",
    "summary(fit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "8d2fa73d-1083-4545-8452-8453e441f3fe",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "            Df Sum Sq Mean Sq F value   Pr(>F)    \n",
       "trt          4 1351.4   337.8   32.43 9.82e-13 ***\n",
       "Residuals   45  468.8    10.4                     \n",
       "---\n",
       "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fit <- aov(response~trt,data=cholesterol)\n",
    "summary(fit)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eff231a2-db82-4697-a62c-afc416d47968",
   "metadata": {},
   "source": [
    "**看着P值是不同的。**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "78ae72e3-ceae-4e52-a69d-d58078641e46",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "            Df Sum Sq Mean Sq F value  Pr(>F)   \n",
       "gesttime     1  134.3  134.30   8.049 0.00597 **\n",
       "dose         3  137.1   45.71   2.739 0.04988 * \n",
       "Residuals   69 1151.3   16.69                   \n",
       "---\n",
       "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fit <- aov(weight ~ gesttime + dose,data=litter)\n",
    "summary(fit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "4157b4d3-3f9d-4a4e-9773-3bcfc9ebe549",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"dataframe\">\n",
       "<caption>A data.frame: 6 × 4</caption>\n",
       "<thead>\n",
       "\t<tr><th></th><th scope=col>dose</th><th scope=col>weight</th><th scope=col>gesttime</th><th scope=col>number</th></tr>\n",
       "\t<tr><th></th><th scope=col>&lt;fct&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;int&gt;</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>1</th><td>0</td><td>28.05</td><td>22.5</td><td>15</td></tr>\n",
       "\t<tr><th scope=row>2</th><td>0</td><td>33.33</td><td>22.5</td><td>14</td></tr>\n",
       "\t<tr><th scope=row>3</th><td>0</td><td>36.37</td><td>22.0</td><td>14</td></tr>\n",
       "\t<tr><th scope=row>4</th><td>0</td><td>35.52</td><td>22.0</td><td>13</td></tr>\n",
       "\t<tr><th scope=row>5</th><td>0</td><td>36.77</td><td>21.5</td><td>15</td></tr>\n",
       "\t<tr><th scope=row>6</th><td>0</td><td>29.60</td><td>23.0</td><td> 5</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A data.frame: 6 × 4\n",
       "\\begin{tabular}{r|llll}\n",
       "  & dose & weight & gesttime & number\\\\\n",
       "  & <fct> & <dbl> & <dbl> & <int>\\\\\n",
       "\\hline\n",
       "\t1 & 0 & 28.05 & 22.5 & 15\\\\\n",
       "\t2 & 0 & 33.33 & 22.5 & 14\\\\\n",
       "\t3 & 0 & 36.37 & 22.0 & 14\\\\\n",
       "\t4 & 0 & 35.52 & 22.0 & 13\\\\\n",
       "\t5 & 0 & 36.77 & 21.5 & 15\\\\\n",
       "\t6 & 0 & 29.60 & 23.0 &  5\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A data.frame: 6 × 4\n",
       "\n",
       "| <!--/--> | dose &lt;fct&gt; | weight &lt;dbl&gt; | gesttime &lt;dbl&gt; | number &lt;int&gt; |\n",
       "|---|---|---|---|---|\n",
       "| 1 | 0 | 28.05 | 22.5 | 15 |\n",
       "| 2 | 0 | 33.33 | 22.5 | 14 |\n",
       "| 3 | 0 | 36.37 | 22.0 | 14 |\n",
       "| 4 | 0 | 35.52 | 22.0 | 13 |\n",
       "| 5 | 0 | 36.77 | 21.5 | 15 |\n",
       "| 6 | 0 | 29.60 | 23.0 |  5 |\n",
       "\n"
      ],
      "text/plain": [
       "  dose weight gesttime number\n",
       "1 0    28.05  22.5     15    \n",
       "2 0    33.33  22.5     14    \n",
       "3 0    36.37  22.0     14    \n",
       "4 0    35.52  22.0     13    \n",
       "5 0    36.77  21.5     15    \n",
       "6 0    29.60  23.0      5    "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "head(litter)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4b5d6c7-52be-4f87-889b-5f4a6ccdba26",
   "metadata": {},
   "source": [
    "**单因素协方差分析：主要用于在控制一个或多个协变量（连续变量，非研究关注的自变量，但与因变量相关）的影响后，检验单一分类自变量（因素）的不同水平对因变量的差异是否具有统计学意义。其核心是通过回归调整，消除协变量对因变量的干扰，更精准地评估因素的效应。**\n",
    "\n",
    "上面的代码中，体重就是因变量，妊娠时间为协变量，药剂为要检验的因素（自变量）。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "036d3805-ed23-4d66-84f6-24f1998b0471",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "\n",
    "协方差分析（ANCOVA）和分层分析是两种不同的统计方法，**二者并不完全等同**，但在控制混杂因素的逻辑上有一定相似性，需结合具体场景区分：\n",
    "\n",
    "\n",
    "**一、核心定义与原理对比**\n",
    "\n",
    "**1. 协方差分析（ANCOVA）**\n",
    "\n",
    "- **本质**：在方差分析（ANOVA）的基础上，加入对**连续协变量**（如年龄、基线值等）的控制，通过回归调整来分离协变量对因变量的影响，从而更准确地评估自变量（分类变量，如处理组、组别）对因变量的效应。  \n",
    "- **模型形式**：  \n",
    "  $$\n",
    "  Y = \\beta_0 + \\beta_1X_{\\text{自变量}} + \\beta_2Z_{\\text{协变量}} + \\epsilon\n",
    "  $$ \n",
    "  其中 $Z$ 是连续协变量（如血压、体重等），通过控制 $Z$ 的影响，分析 $X$ 对 $Y$ 的净效应。\n",
    "\n",
    "**2. 分层分析**\n",
    "\n",
    "- **本质**：将数据按某一**分类变量**（如性别、地区、疾病严重程度等）分成不同亚组，在每个亚组内分别进行分析（如分层后的卡方检验、分层回归等），以考察自变量在不同亚组中的效应是否一致，或控制分层变量的混杂作用。  \n",
    "- **核心目的**：  \n",
    "  - 评估自变量效应的**亚组异质性**（如处理效应在男性和女性中是否不同）；  \n",
    "  - 控制分层变量的混杂影响（如按年龄分层后，分析吸烟与肺癌的关系）。\n",
    "\n",
    "\n",
    "**二、关键区别**\n",
    "\n",
    "| **特征**         | **协方差分析（ANCOVA）**                          | **分层分析**                                  |\n",
    "|------------------|---------------------------------------------------|-----------------------------------------------|\n",
    "| **控制变量类型** | 连续变量（协变量，如身高、基线血压）              | 分类变量（分层因素，如性别、组别、地区）      |\n",
    "| **模型形式**     | 单一模型中同时纳入自变量和连续协变量              | 分亚组建模，或在模型中加入分层变量与自变量的交互项 |\n",
    "| **核心目标**     | 调整连续混杂因素的影响，估计自变量的“净效应”      | 考察自变量效应在不同亚组中的差异或控制分类混杂因素 |\n",
    "| **应用场景**     | 例如：比较不同药物对患者疗效时，控制基线血糖（连续变量）的影响 | 例如：分析教育水平对收入的影响时，按性别（分类变量）分层分析 |"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9bdcaaa3-6f60-4c40-8d43-554d6bf15344",
   "metadata": {},
   "source": [
    "下面以一个具体的例子来详细说明单因素协方差分析的计算步骤，以及 F 值和 P 值的计算方法。\n",
    "\n",
    "**案例背景**\n",
    "假设我们要研究三种不同的教学方法（A、B、C）对学生数学成绩（因变量 $Y$）的影响，同时考虑学生的入学数学成绩（协变量 $X$）可能会对最终成绩产生影响。我们收集了 30 名学生的数据，每种教学方法各有 10 名学生，具体数据如下表所示：\n",
    "\n",
    "| 教学方法 | 入学成绩 $X$ | 最终成绩 $Y$ |\n",
    "| --- | --- | --- |\n",
    "| A | 70, 72, 75, 78, 80, 82, 85, 88, 90, 92 | 75, 78, 80, 82, 85, 88, 90, 92, 95, 98 |\n",
    "| B | 72, 75, 78, 80, 82, 85, 88, 90, 92, 95 | 78, 80, 82, 85, 88, 90, 92, 95, 98, 100 |\n",
    "| C | 75, 78, 80, 82, 85, 88, 90, 92, 95, 98 | 80, 82, 85, 88, 90, 92, 95, 98, 100, 102 |\n",
    "\n",
    "**单因素协方差分析的计算步骤**\n",
    "\n",
    "**步骤 1：计算各项均值**\n",
    "\n",
    "- **总均值**：分别计算协变量 $X$ 和因变量 $Y$ 的总均值 $\\bar{X}$ 和 $\\bar{Y}$。\n",
    "- **组内均值**：计算每个组（教学方法）内协变量 $X$ 和因变量 $Y$ 的均值 $\\bar{X}_i$ 和 $\\bar{Y}_i$（$i = 1,2,3$ 分别代表 A、B、C 组）。\n",
    "\n",
    "假设经过计算得到：\n",
    "$\\bar{X}=82$，$\\bar{Y}=88$\n",
    "\n",
    "$\\bar{X}_1 = 80$，$\\bar{Y}_1 = 86$\n",
    "\n",
    "$\\bar{X}_2 = 83$，$\\bar{Y}_2 = 89$\n",
    "\n",
    "$\\bar{X}_3 = 86$，$\\bar{Y}_3 = 92$\n",
    "\n",
    "**步骤 2：计算回归系数 $b$**\n",
    "\n",
    "协方差分析假设协变量 $X$ 和因变量 $Y$ 之间存在线性关系，通过最小二乘法计算回归系数 $b$。计算公式为：\n",
    "\n",
    "$$b=\\frac{\\sum_{i = 1}^{k}\\sum_{j = 1}^{n_i}(X_{ij}-\\bar{X}_i)(Y_{ij}-\\bar{Y}_i)}{\\sum_{i = 1}^{k}\\sum_{j = 1}^{n_i}(X_{ij}-\\bar{X}_i)^2}$$\n",
    "\n",
    "其中 $k$ 是组数（这里 $k = 3$），$n_i$ 是第 $i$ 组的样本量（这里 $n_i = 10$）。\n",
    "\n",
    "假设计算得到 $b = 0.8$。\n",
    "\n",
    "**步骤 3：计算调整后的因变量均值 $\\bar{Y}_{i_{adj}}$**\n",
    "\n",
    "调整后的因变量均值是在控制协变量影响后的因变量均值，计算公式为：\n",
    "\n",
    "$\\bar{Y}_{i_{adj}}=\\bar{Y}_i - b(\\bar{X}_i - \\bar{X})$\n",
    "\n",
    "对于 A 组：$\\bar{Y}_{1_{adj}}=86 - 0.8\\times(80 - 82)=87.6$\n",
    "\n",
    "对于 B 组：$\\bar{Y}_{2_{adj}}=89 - 0.8\\times(83 - 82)=88.2$\n",
    "\n",
    "对于 C 组：$\\bar{Y}_{3_{adj}}=92 - 0.8\\times(86 - 82)=88.8$\n",
    "\n",
    "**步骤 4：计算平方和**\n",
    "\n",
    "- **总平方和 $SST$**：反映因变量 $Y$ 的总变异，计算公式为：\n",
    "\n",
    "$$SST=\\sum_{i = 1}^{k}\\sum_{j = 1}^{n_i}(Y_{ij}-\\bar{Y})^2$$\n",
    "\n",
    "- **协变量平方和 $SS_{X}$**：反映协变量 $X$ 对因变量 $Y$ 的影响，计算公式为：\n",
    "\n",
    "$$SS_{X}=b^2\\sum_{i = 1}^{k}\\sum_{j = 1}^{n_i}(X_{ij}-\\bar{X})^2$$\n",
    "\n",
    "- **组间平方和 $SSA$**：反映自变量（教学方法）在控制协变量后对因变量的影响，计算公式为：\n",
    "\n",
    "$$SSA=\\sum_{i = 1}^{k}n_i(\\bar{Y}_{i_{adj}}-\\bar{Y})^2$$\n",
    "\n",
    "- **残差平方和 $SSE$**：反映随机误差和未被模型解释的变异，计算公式为：\n",
    "\n",
    "$$SSE = SST - SS_{X}-SSA$$\n",
    "\n",
    "假设经过计算得到：$SST = 500$，$SS_{X}=100$，$SSA = 80$，$SSE = 320$\n",
    "\n",
    "**步骤 5：计算自由度**\n",
    "\n",
    "- **协变量自由度 $df_{X}=1$**\n",
    "- **组间自由度 $df_A=k - 1=3 - 1 = 2$**\n",
    "- **残差自由度 $df_E=N - k - 1=30 - 3 - 1 = 26$**(本质上是总样本量减去模型中估计的所有参数（组截距 / 组效应 + 协变量斜率）的结果。)\n",
    "- **总自由度 $df_T=N - 1=30 - 1 = 29$**\n",
    "\n",
    "**步骤 6：计算均方**\n",
    "\n",
    "- **协变量均方 $MS_{X}=\\frac{SS_{X}}{df_{X}}$**\n",
    "- **组间均方 $MSA=\\frac{SSA}{df_A}$**\n",
    "- **残差均方 $MSE=\\frac{SSE}{df_E}$**\n",
    "\n",
    "计算可得：$MS_{X}=\\frac{100}{1}=100$，$MSA=\\frac{80}{2}=40$，$MSE=\\frac{320}{26}\\approx12.31$\n",
    "\n",
    "**F 值和 P 值的计算**\n",
    "\n",
    "**F 值的计算**\n",
    "\n",
    "F 值用于检验自变量（教学方法）在控制协变量后对因变量的影响是否显著，计算公式为：\n",
    "\n",
    "$$F=\\frac{MSA}{MSE}$$\n",
    "\n",
    "将上面计算得到的 $MSA = 40$ 和 $MSE\\approx12.31$ 代入公式，可得：\n",
    "\n",
    "$F=\\frac{40}{12.31}\\approx3.25$\n",
    "\n",
    "**P 值的计算**\n",
    "\n",
    "P 值是在原假设成立的情况下，得到当前 F 值或更极端值的概率。可以通过 F 分布表或统计软件来计算 P 值。\n",
    "在 R 语言中，可以使用 `1 - pf(F, df_A, df_E)` 来计算 P 值，其中 `F` 是计算得到的 F 值，`df_A` 是组间自由度，`df_E` 是残差自由度。\n",
    "\n",
    "```R\n",
    "F <- 3.25\n",
    "df_A <- 2\n",
    "df_E <- 26\n",
    "p_value <- 1 - pf(F, df_A, df_E)\n",
    "p_value\n",
    "```\n",
    "\n",
    "假设计算得到 $P\\approx0.053$\n",
    "\n",
    "### 结果解释\n",
    "- **F 值**：F 值越大，说明自变量（教学方法）在控制协变量后对因变量（数学成绩）的影响越显著。\n",
    "- **P 值**：通常以 $P = 0.05$ 作为显著性水平的临界值。如果 $P < 0.05$，则拒绝原假设，认为自变量在控制协变量后对因变量有显著影响；如果 $P\\geq0.05$，则不能拒绝原假设，认为自变量在控制协变量后对因变量的影响不显著。在这个例子中，$P\\approx0.053\\geq0.05$，说明在控制入学成绩的影响后，三种教学方法对学生数学成绩的影响不显著。 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "5d07cc06-2afa-4879-a3f8-a55a6202f5a0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "            Df Sum Sq Mean Sq F value  Pr(>F)   \n",
       "dose         3  109.9   36.64   2.196 0.09628 . \n",
       "gesttime     1  161.5  161.49   9.679 0.00271 **\n",
       "Residuals   69 1151.3   16.69                   \n",
       "---\n",
       "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fit <- aov(weight ~ dose+gesttime,data=litter)\n",
    "summary(fit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "c46f2afd-5ba7-4cfb-8255-8304ac98d055",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "            Df Sum Sq Mean Sq F value   Pr(>F)    \n",
       "supp         1  205.4   205.4  12.317 0.000894 ***\n",
       "dose         1 2224.3  2224.3 133.415  < 2e-16 ***\n",
       "supp:dose    1   88.9    88.9   5.333 0.024631 *  \n",
       "Residuals   56  933.6    16.7                     \n",
       "---\n",
       "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fit<-aov(len~supp*dose,data=ToothGrowth)\n",
    "summary(fit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "8ade6d4d-0446-408d-a488-50bfdb9f34a9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "60"
      ],
      "text/latex": [
       "60"
      ],
      "text/markdown": [
       "60"
      ],
      "text/plain": [
       "[1] 60"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nrow(ToothGrowth)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "c7b69502-b225-4063-9fe5-b8947f3e0053",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "package 'boot' successfully unpacked and MD5 sums checked\n",
      "\n",
      "The downloaded binary packages are in\n",
      "\tC:\\Users\\xie.xiaokang\\AppData\\Local\\Temp\\RtmpYRYtSj\\downloaded_packages\n"
     ]
    }
   ],
   "source": [
    "# install.packages(\"boot\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "9269ca6e-904d-4abb-9c16-0730a95c5de1",
   "metadata": {},
   "outputs": [],
   "source": [
    "rsq <- function(formula,data,indices){\n",
    "    d <- data[indices,]\n",
    "    fit <- lm(formula,data=d)\n",
    "    return (summary(fit)$r.square)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "3c09a620-bcf4-43b1-9696-680b77b4c9cc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "ORDINARY NONPARAMETRIC BOOTSTRAP\n",
      "\n",
      "\n",
      "Call:\n",
      "boot(data = mtcars, statistic = rsq, R = 1000, formula = mpg ~ \n",
      "    wt + disp)\n",
      "\n",
      "\n",
      "Bootstrap Statistics :\n",
      "     original     bias    std. error\n",
      "t1* 0.7809306 0.01379126  0.05113904\n"
     ]
    }
   ],
   "source": [
    "library(boot)\n",
    "set.seed(1234)\n",
    "results <- boot(data=mtcars,statistic=rsq,R=1000,formula=mpg~wt+disp)\n",
    "print(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3943e3e9-005d-4271-ae1f-71c983e42e5e",
   "metadata": {},
   "source": [
    "- `original`：\n",
    "    - 其值为0.7809306，指的是基于原始数据集计算得到的统计量的值。在此处，它是依据`mpg ~ wt + disp`这个线性回归模型算出的决定系数（\\(R^{2}\\)）。\\(R^{2}\\)反映了回归模型对数据的拟合优度，取值范围在0到1之间，越接近1，说明模型对数据的拟合效果越好。\n",
    "- `bias`：\n",
    "    - 其值为0.01379126，代表自助法估计的偏差。偏差的计算方式是自助法样本统计量的均值减去原始样本统计量的值。如果偏差值接近0，表明自助法估计的统计量和原始样本统计量较为接近，没有明显的偏差；若偏差值较大，则说明自助法估计可能存在偏差。\n",
    "- `std. error`：\n",
    "    - 其值为0.05113904，指的是自助法估计的标准误差。标准误差衡量的是自助法样本统计量的变异性，也就是不同自助法样本计算出的统计量之间的波动程度。标准误差越小，表明自助法估计越稳定、可靠。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "2cc67488-38bb-41eb-9df5-d6c914de16ee",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAANlBMVEX9/v0AAABMTUxnaGd7\ne3uLjIuZmpmmpqaxsrG7vLvFxsXOz87T09PX2Nff4N/n6Ofu7+79/v03sFGxAAAAEnRSTlP/\n/////////////////////wDiv78SAAAACXBIWXMAABJ0AAASdAHeZh94AAAgAElEQVR4nO2d\njYKqqhpAD5ZZ0670/V/2+FtQIoioCGvdc/c0acBHrEEQ9b8KABbz394FAIgBRALwACIBeACR\nADyASAAeQCQADyASgAcQCcADiATgAUQC8AAiAXgAkQA8gEgAHkAkAA8gEoAHEAnAA4gE4AFE\nAvAAIgF4AJEAPIBIAB5AJAAPIBKABxAJwAOIBOABRALwACIBeACRADyASAAeQCQADyASgAcQ\nCcADiATgAUQC8AAiAXgAkQA8gEgAHkAkAA8gEoAHEAnAA4gE4IFwRRJCqK8+b8hcNinMNRNC\nzmmbXOE4HFykRzYml3eudd6SSBvlCgfi4CKN91LeOQnx3D5XOBAHEml6p43KsmWucCAOJFL/\no7ye61f5X/eOGHa7X5qjr3v/kVf92/kmffJ1EkX96i+vX5+K15De7SROj6q6ZeL8ULNX0lPM\nkXMF6DicSK+sb8hnpUmf+9d5+4lHv8vnk6f2A8Ne4tG92/3+Kt7vvZHT+zIHkeCXw4lUdxR1\nZ1TWDf0mNel8MKQzKXv/OnxSNB+71eqUVVUofgiRyQ52KOkhEhgJWqQfHbp/mwOzsu5k3u9V\n9/rnrayP+uqf9dHYXy1H8yP7fLIRqJkzeCkp1e/emu7q2f745P2VHmMkMHE4kRo53kOhoUlf\nmv6poWhnqfOu+TdCDZ+8fyXd/ftQfnx2+EoPkcDE4US6dm/0Ln02le3vr/aNbGjo35vrHf6K\ns3iLVP38eH9OTg+RwETQIqmv+h/FMLJ5/WwaXolfkbrf/06SmdMiKa8QCQwcT6Sq/Oum1M7V\naI+UjfZI7a/Nod7pcnvO6pGy742IBL8cUKSG9izP573cOEZqt576940i5YyRYB6HE+nUdxaf\nrqLUztqJL0n6n+YeyWLWrqwAPhxOpLqNn1/tnEOzUqGZw2t+vs+0dktLf88jtQmd253vmVGk\n7/S+RBpyBRg4nEjvyYZmiNTMU3cvznK7b3sUoaxsaN9+DPMU7Yz3lEhf6X1tfOcK0HM8kbrx\n0bkbwzSjmd6dSyadYHo2a+3uP5I0b2eX52tYsDCS+oCS3vfGXL06CSBckTxQdgMpgNWJUiTR\nrsernmd1AR3AakQp0meqQF0aBLAWUYr0vtSCuTXYiChFqsprcx1EdqE/go2IUySAjUEkAA8g\nEoAHEAnAA4gE4AFEAvAAIgF4AJEAPIBIAB5AJAAPIBKABxAJwAOIBOABRALwACIBeACRADyA\nSAAeQCQADyASgAcQCcADiATgAUQC8AAiAXgAkQA8gEgAHkAkAA8gEoAHEAnAA4gE4AFEAvAA\nIgF4AJEAPIBIAB5AJAAPpCWS2LsAoSJAYX4NIhJU6VaMJm5EMpBqezFCxSggkgHaiwYqRgGR\nwA1EUkAkcCNRkXRhIxK4kaZI2qgRyUCa7cWCNCsGkVxJs71YkGTF6INGJANJthcbkqwYRHIm\nyfZiAxWjgEgGaC8aqBgFRAI3EEkBkcCNBEWaChmRwA1Est6mIS2REmwvdqRXMZMRI5KB9NqL\nJelVDCItIb32Ykl6FYNIS0ivvVhCxSggkgHaiwYqpr3c/v1y/qfTEgk0IFJr0aASIoEbqYn0\nG6+Q30ckcCMxkXQeffk0h7REcr7ZUuwkViWItBDxryGxVmNDWlUyEi0izQKRNKRVJWPRMkaa\nAyJpoEqYtZsDImmgSjiPNAdE0kCVKCCSAUTSkFSVmINFJAOIpCGlKrGIFZEMIJKGlKoEkZaD\nSBoSqhKbUBHJACJpSKhKEMkDiKSBKlFAJAOIpIEqUUAkA4ikgSpRQCQDiKQhmSqxCxSRDCCS\nhmSqBJG8gEgaUqkSyzgRyQAiaUilShDJD4ikIZUqQSQ/IJIGqkQBkQwgkgaqRAGRDCCSBqpE\nAZEMIJKGNKrEOkpEMoBIGpKokp8g24vLx27PhkgGEElDElXyFeTndic/KiGSAUTSkEKVSDF+\n7hPa/6vf1RZEgioxkWp5Wo+a//q3vzur+akjElRpiPRGfP5DJDcQSUNKVdLL0000IJITiKQh\npSqR5BGMkdxITaTHNW+H1XnxmN4xgSqRb6M6zDMIZu3cSEuk8iQ+nCd3jb9KhPKy/4/zSI6k\nJVIhsr9n++p1z0QxtWv8VaJMfr9n7cw724JI8ZKJ5/v1U2RTu0ZfJWqApufNIZKBtERSmsr0\ncwqjrxK1KkzhIpKBtESiRxpjfHrha5/5ySJSvNRjpPurfcUY6YOQ/p3eZxaIFDFnadbuVE7t\nmUyVvI/rfB/pIlLMPIr2PFKWX9M+jyQtsxPDwZ3V/vYgElSxiySkV4NEiLQERNIQd5XIixr6\n/zNGWkRqIj2Lbph0yv+md4y6SoTySreaYfwT1iBSxFylyYZ8cs+oq0R8veI80mLSEukuLq+q\nepzz6nk7ifvUrklUyXs9gzFaRDKQlkhn0U55P8W11mm6S0qgSqTeGZGWkpZIQ4NpFzUkvkRo\nuPbIRiNEMpKWSFnXI5UWi2JirpLhBOxwDZ+FSYhkIC2RCnF+VNUrF5eqvNT/TBBxlQzdcjWM\njxBpOWmJNCwRysq68WSv761izqjhuAwive8bhEjLSUyk6lardLrWL7JicqldxD3SEFmvkZVH\niGQiNZGsibdKpA7JvvdFJAO9SNEfzswm9rro51tsv3VEMtCL9I+O6YvY6+J91y27QLcU6XUR\nWX30fTuJbPKSsaBIWKR0p7+rj0G2RyEbilRmTTd5u1rc6ikgEEm3datibMz33Hd4IhXNpctF\nJi5lVRbTlzEHRMIiTRNrXbynGirbRQ3yx+bgKlLWd5XtrOr0jTUCApE0RFoX0tx3Zd8fbSqS\nUrTfEgZ6sg+RNERaF9KRXRWoSJkkUnmYWz2lJlLq9/4Wnx9CmC6L/f3cHJaOkZoz5oYxUkBf\nUloice/vns86O+sPzGaDWbuAvqS0ROLe3z3DXIP9B+bnscF5pIC+pLRE4k6rPfPmvqtQVzYE\n9CWlJVLq9/4W0s9ZQyREMpGWSIn3SOL7BSJ5Iy2REr/3tyTSzLkGRDKRlkhp3/v7E9LsuQZE\nMpGYSEnf+1sW6X3H79mftQaRoArqO/JNf435nLkGRDKBSBrirYvBI0TyCSJpiLYuPv0Rh3Ye\nQSQN8dXFsJy674zmrZ1GJAOIpCG6uujXVA/XH8y8CAGRDCCShujq4i3S26bZn54HIkEV1Hfk\nBbVDQiTfIJKG2OpiEKlyOK6rEMkIImmIsi7cL9JGJAOIpCHGulhwswNEMoBIGmKsC1Eh0log\nkoa46uJrqsExgXkgElRBfUceeE81OMwzfBKYByJBFdR3tJx+UcPwqEvnFGaBSFAF9R0tp+uQ\nKucDO0QygkgaoqoLdVWDcwrzQCSogvqOPNH3SG6RIZIBRNIQXV2816q6fXr+RxAJqqC+Iy8s\nmfuuEMkIImmIqC76O52IBUd2iGQCkTTEUxff10+4JzIPRIIqqO9oIUqHhEgrgUgaoqkL+foJ\nZ48QyQQiaYimLrx0SIhkApE0RFUXwzyDs0eIZAKRNERVF0v7I0QygkgaYqqL5R4hkglE0hBL\nXagnkRYlMxNEGtlr8RdxOCKJVZ1qWJbOTBBJt1dSPVYksQ4iVYi0LoikIY5YfXVIiGQCkTTE\nEask0qzHIY0mNBNE0u2FSEdE+DiyQyQTiKQhmli9eIRIJhBJQzSx9g4h0rogkoYYYvU21YBI\nRrQiKdWPSIdEmftentRMEEna8A+RjovHDgmRTCCShghiFR4un/ikNRtEkjYg0oGRns2HSGuD\nSBpiiNVfh4RIJhBJQwyxItJ2IJKGCGL16BEimTCK1IFIh0P49AiRTBhF0m2PnaPH6tcjRDKB\nSBqOHqtfjxDJBCJpOHisnjskRDKBSBoOHqtfjRDJCCJpOHisiLQxiKTh2LH69giRTCCShkPH\n6nmAVCGSEUTScORY/XuESCYQScOBY13BI0QygUgaDhwrIu0AImk4bqxreIRIJhBJw3FjRaQ9\nQCQNh411FY8QyUSSIt1OQuT36X0OG+sqHiGSlqGykxKpa1vnLvJietdNCrQCq3iESFpMosQr\nUiGKsqpehbhN7rpRkXyzjkeIpCVdkTJRNq9LcZrcdZMC+QeRNiZdkYYmNt3UDhrrSh4hkpZ0\nRboMImWTu25QnBVApK1JVKT8eruLv/plWUzPNhwz1nWm7CpE0pOoSO9GJkRWTu66UZG8sppH\niKRlVBTpO4hSpOr5vN3yvJ1yKCY9CuI7mg0ibc+oKJIucYpkzxFjXc8jRNKCSNMcMVZE2oE0\nRXpdRHZt1whl0wsbgviO5rKeR4ikJUmRyqxpZbdr29jOk7seMNYVPUIkLUmK1E55F5m4lDFO\nfyPSHiQpUtZGILo1QrGdkF3TI0TSkqRIQnz+jW6JECLtQpIiZZJIZVw9kjRjh0hbkqRIwxip\nORkb2RgJkXYiSZEMs3Zi5da4ImuXHJF0JClStOeRVv8TgEg60hTJnmPFunpXikg6EGmaQ8W6\n/iEpIulApGkOFSsi7QcixXMeaYMpEkTSgUiINCeL+R9BJP32ncu8JQeKdYs5e0TSgUjTHChW\nRNoTRJrmOLFuchIZkXQkKtLjmrdNLi8e0zseJ1ZE2pUkRSpPUquL5cK+TVY1IZKOJEUqRPb3\nbF+97lkki1Y38QiRtCQpUiae79fPSC6jQKR9SVIkpbFFch4JkfYlSZEi7JG28QiRtCQpUj1G\nur/aV9GMkRBpZ5IUaXhYX8spint/b+MRImlJU6TqUbTnkbL8Gsd5pI06JETSkqhI1hwj1k1O\nxrYZzf9I9CL1FY9Ikxwi1s08QqSx3EfcQaRvDhErIu0JItlwhFi38wiRxnJHJAsOEKtApPWz\nmModkSw4QKwbeoRIY7kjkgXhx7plh4RIY7kjkgXhx7qlR4g0lrtJJP30+K7l3pbgY93UI0Qa\ny90kkn77ruXeluBjRaS9vyREsiH0WLf1CJHGckckC0KPdVuPEGksd0SyIPBYN+6QEGksd0Sy\nIPBYN/YIkcZyRyQLAo8VkRoQKXwCj3VjjxBpLHdEsiDsWLfukBBpLHdEsiDsWBGpBZHCJ+hY\nt56zQ6TR3BHJgqBj3dwjRBrLHZEsCDnW7TukTUUqi+aOg9eTEOe/6T0RKXwCjnUHj7YU6ZXV\ncZVZF2DITzpAJBsCjjVykS4iL+t/Lq/aqUvId/FcINLG396ehBvlHh5tKZIQZf9PfZQX8n2l\nF4iUUL8UbpTxi1Q1d2mXfvnavEsFjIBINgQb5T7NaNNDu2dVXbvHHZTTgyRECp9go4xepKfI\nimeVZ7VJ95O4ey6VRxDJhlCj3Om4Zsvp73v2ifE6uScihU+oUSYgUlX9XdqH/ebX1/R+iBQ+\ngUa510CblQ0juSOSBYFGiUgyiBQ+gUa518wvIo3kjkgWhBnlbmdQEGkkd0SyIMwoEUkBkcIn\nyCj3O6WPSCO5I5IFQUaJSCqIFD4hRrnjIjNEGskdkSwIMUpE+gKRwifAKHf0CJHGckckCwKM\nEpG+QaTwCTDKHT1CpLHcEcmC8KLc0yNEGssdkSwIL0pE+gGRwie8KBHpB0QKn+Ci3NUjRBrL\nHZEsCC5KRPoFkcIntCj3nPuuEGk0d0SyILQo9/UIkcZyRyQLAouy9geRfkCk8Aksyp07JEQa\nyx2RLAgsyp09QqSx3BHJgsCiRKQx9uqeexDJgrCi3PvIDpGUbA3uIJJEWFEi0iiIFD5BRbnz\nSaQKkdRsEcmeoKLc3SNEUrJFJHuCihKRxkGk8Akpyv09QiQlW0SqHte8bZF58ZjeMaQoEUkD\nIu1EeZKG7SE/MFtlf48QSck2eZEKkf21j1SsXvcs5AdmKwTQISGSkm3yImXdk0lbniE/MFsB\nkXQg0k4oTXG6XQYUJSLpQKSdOHiPtGcZ5n8EkeIVqR4j3bvHkh5ojBRCh4RISrbJi1SdpVm7\nUzm1ZzhRIpIWRNqNR9GeR8ry61HOI+2/zq4txfyPIFLMIlkTTJRBeIRISraIZE8oUYbRISGS\nki0ivS4iu1bV7SSyyamGEEXatxjzP4JI8YpUZk2LvF2Ps0QokA4JkZRskxepaKa8i0xcyqos\nDjH9HYhHiKRkm7xIWVt2IdqJ70OckEWkKRBpJ4T4/HuMJUKBeIRISrbJi5RJIpVH6JFC6ZAQ\nSck2eZGGMVJR9q/1hBElIk2CSDthmLUTMjsU74dwioNIcrbJi3S080jBeIRISraIZE8QUSLS\nNIgUPiFEGY5HiKRki0j2hBAlIhlApBAI/zxSOB4hkpItIskEL1JAHRIiKdkikj0BRIlIJhAp\nfPaPMpyTSBUiqdkikj37RxmSR4ikZItIR7r3NyIZQaSdONK9v4PyCJGUbJMX6UD3/g5qhIRI\narbJi3SgO60ikgWItBPHufd3YB4hkpJt8iIdp0eKTKTT9eWtKBKItBOHufd3aB4tFamOYw2X\nEGkvjnLv79A8WipS+XdZwyVE2o2D3Ps7NpEaHteTb5cQKXx2jTI4jzxNNjybS/1vy0szgEjh\ns2eUwY2QPIl0P1ucCZ8DIoUPIiksF6m81t3R6V7WNuV+yoRIR2DHKAP0aLFIj2ayoehOPvgL\nC5HCB5EUlp5Hqjuj2zBLOn0Cbw6IFD5BiLRfGb5Zeh4pv3sripzsGolaZItI9uwXZYgd0uLz\nSN4KooBI4ROCSLsV4ZfFKxv6F5m3w7o2WZ+JzcgWkezZLcogPfIl0stvVIgUPoiksECku5A5\n7VwqL9kikj17RRnkCGlZjyRflnwyrMxavVReskUke/YXaacCjONrjOQXRAqfvb6jOEVaCUQK\nn91F2id/HQtEakJZ6e8DIoXPPlGG2iEhkpItItmzS5TBesShnZItItmDSAqIJGeLSPbsEWW4\nHi0W6XaqqtfJ8+w3Ih0ARFJYKNK9iah9EjbnkRBp9SzD9WipSGfxVz3FqfrzeHlshUhHAJEU\nPJyQfTb3P2PWDpFWzzFgj3yIlIs7IiHSBjkG7NHyQ7vnvbkwlkM7RFo9w5hFaleAX5sgvV4p\ni0i+eJ1XuvZyT5G2ztmGxdPf3R2iT3+eytOBSL4oxeN+XSVlRFLghKycbXwiVflFXFZJeD+R\nts7YCkSSs41OJCFOPu+AqyS9TrITGSLSbBDJG/kli6NHCvvIbrFI19Ma4SGSL0rxjGSMFLdI\n13XCi0KkIL72l9fTEjK7ibRxvpYsFClb5wg8DpGC6568snFggXdI3LNByTZKkeK492CrULwi\n5WKV832I5JEo7j0Y9rKGarFIr+zs90qkDkTyRDT3HgzdIw8PY2ayIWCRorn3ICI5gUgeiWEc\nG7xHnJBVso1SpJVAJAVEkrNFJHsQSWGxSPe8vbjv5ak8HYgUPlsGFr5Hi0U6d9GJzKtJiBQ+\niKSwUKSbOJdNdDe/a/URKXw2DOwAHi1fIlR2s0LM2iHSelnFL1J3/+8KkRBpxZwSEOnU90jP\nI581/2SLSPbsIdJmWc7Hzxjp7nkVOCKFz2aBHaJDWjxrl/ch+r3sBZHCZweRtsrRBS/nkUTu\n9yZCiHQAtgrsGB0SKxuUbBHJnu1F2ihDNxBJzhaR7EEkhWUi3S/NMv1z4fuaJEQKn81F2ig/\nR5aI9Dq/gzyz1g6R1skmfpHKTJzuzZXmr7+T8Ho7AEQ6AIiksECkQprzPguvN09DpPDZJrCj\neLREpJP4HM+9eKwLIq2SSwIiKbHNCtS4MyKFzyaBHeQkUoVIaraIZM/GIm2R2xI2FEmoeC6V\nDxBpBlsEdpwOaUuRHhkiHYTbSYjc8AzGbUXaILNlLBLJWoyOMu/PN43uPDOxNUCkrub784PF\n9K5blAaRxvkT4q+Kf4zUc1SRClGUVfUqpi+N2VSk9fNaysZr7V5nkZfRizTSPR2E9pvJuhu6\nl9MXayKSwuaLVq8iuyNSqCi3Ddh7HHskj3ZY/f08mWsGkXai/WYug0iT674QSWGPyyguiBQq\nQuTX270dyJbF9GwDIilwPZKcLSJ9Wq4Q2eSzr1YP7FAeIZKSbfIiVc/n7Zbn7ZRDMf0MOURS\nSF4k+btCpBlsJ9LaGXkBkaTmjkgzQCQFREIkCZFZ3zUAkRQQCZEk6mabWz5ee+3AjjVEQiRE\nkhHingnDLMOw6+pF6SxCpAUg0k7UrbbMhbgYVn63u65dkkN1SIiESDJts322N8+9PXed/kYk\nHyDSTvTN9llkxja8lUjrZuMNREIkiU+7fd7yEyLZg0iIJGFot0Jm5YIgkgcQaSdmtFtEUkAk\nRHJj1cAO5xEiIZIjiKSASIjkxjYirZmJVxAJkXTsOP2NSH7YXqQORJLZT6TjHdkhkosbiYg0\nDSIpIBIiubFiYFudrPIJIiGSG5uItF4evkEkRFJ4XPO2CeemBwOvF9gROyREQiSZ8iS14uln\nx20h0mpZ+AeREEmiENnfs331aq7wm9oVkRTSFWn4shBJIhPP9+vnTndaPaRHKYvk7ka8Is14\n5hUiKSASIkkE0CMdcqoBkRBJoR4j3btn1e82RjqmR4iESApnqUM47XLvb0TyCCLtxqNozyNl\n+XWf80gH9QiREMkRRFJAJERyA5EUEAmR3FhbpHWSXw1EQiQ31gnsqB0SIiGSI4ikgEiI5AYi\nKSASIrmxskirpL4iiIRIbqwS2GE7JERCJEcQSQGREMmNNQI76ILVBkRCJDfWFWmFxNcFkRDJ\njRUCO3CHhEiI5MiqIvlPe20QCZHc8B9YpxAieQSRwmclkQ7qESIhkiOIpIBIiOSG98AO7REi\nIZIjiKSASIjkxnoi+U54ExAJkdzwHdhnzs5zwtuASIjkxioi0SP5BZHCx3NgR17V0IBIiOTG\nCiId9WRsAyIhkht+AxMHHyIhEiI54luk6sgHdoiESK6s0SMhkmcQKXxWEslrqhuCSIjkBiIp\nIBIiueE1sMMf2SESIjniM7Cjn0SqEAmRXPEY2Hvi+7geJSlS/7cPkRbhV6Sjd0hpirTUDUSq\nfH5HApHWApHCx6dI1eE9QiREcmSVHslbmpuDSIjkhl+RDrzKrgOREMmNFXokbynuACIhkhuI\npIBIiOQGh3YKiIRIbngL7N0fHbqqEAmR3PApUjcF7ivBXUAkRHLDV2D96OjgR3aIhEiOeBSp\nOvLNGnoQCZHc8CZSt2b14Ed2iIRIjngKTEQxZ4dIiOSKR5Ei8AiREMkRXyId/O5BA4iESG74\nCSyG9aotiIRIbvgTqeLQbi0QKXy8iRTB3HeFSIjkis8e6cg3a+hBJERyw/MYyUtqO4JIiOQG\nIikgEiK54fHQLgKPEAmRHPESmGCyYVUQKXw8iTRM3PlIbU8QCZHc8BHYMDyiR1oJRNqNxzVv\nW3dePKZ39CdSDEd2iIRIMuVJusnceXJXD4FFcUO7HkRCJIlCZH/P9tXrnolialevIi1Pa28Q\nCZEkMvF8v36KbGpXTyLFMWeHSIikoDTp6fa9PLDPRebHqqRREAmRJDbtkYabcMXgESIhkkw9\nRrq/2lfrj5FimmpAJERSOUvN+1RO7YlICoiESAqPoj2PlOXXlc8jRfFUpA+IhEhuLBfpPWnn\npTw7g0iI5IYPkWKZ+64QCZEURGY4oJN2XZyXiGbuu0IkRFJoFtlNTjFIuy7PKpL1qi2IhEgS\nQjSz3lYq+RkjRXIWCZEQSaFu1WUuxOVuseuyjCJ5mssbREIkibZ7eDYT4PntOd0xLQpMxHZk\nh0iIJNM362eRGacBlopURTT3XSESIil82vXzlp9WE0l0HxfxHNkhEiLJGDoIIbMkm2HyO5ap\nBkRCJIUZ7drTCdmFqQQDIiGSGwsDi2mdXQMiIZIbfkRijLQqiBQ+CwIT8mjLX4l2BZEQScda\n099iSDumQRIiIZKOFUXqj+ni6ZAQCZEccQ6sNyiS5yINIBIiubFMpEgewfwBkRDJjWVLhCKb\nakAkRPpim3t/R+cRIiGSzEb3/o6vQ0IkRJJZ/97f7cfet1h1SyNEEAmRJNa/06o0V4dIjpSX\n+nChv/Zy7ftKT6aOSDpWv/d33yENJjmlESQbilR2F4vlXb7JiHSo0cDqPdIgUmSnYzcVqRC3\n2qZb1o5h0xHp35E6ptXv/f3ukQ7018WKDUXKurxe2emFSKGyzb2/+2uRjlIpNmwo0lBv5fk8\nJpKvqy8tCoJIera49/exDnft2FCkkxj+wp3O9EiHZ/llFP7KEgAbinQTl/7VS5wR6ei4RCSU\nT0ZVJ1tOfxdve+6Gv0eIFD4OEQn1k1HVyaYnZJ/58Op1QaSDs0ik9g9pVHXCygZEcmN+RO9P\nRHYpUgsiIZIbiKSASIjkhntEHNp1IBIiVQtFWphCgCASIrmxQKRhrZ2/wuwPIiGSG7MjendE\nw9VIngu0L4iESG7MjUiaaojqwUg9iYikLO5CJB+4ijRcQYFIhxRJbs2I5IOZEalz39LPSEAk\nRHLDXaTh2M53iXYFkRDJDceIort7fg8iIZIbbhFtdc3Z5iASIrnhLNLn6C4mEAmR3JgVkXQO\nScg/4wGREMkNZ5GiPB+LSIjkyJyI1Cm7CM8iIRIiueIkkmSS/xLtCiIhkhvuPVKEGiESIrni\nFNEwPKI6KkRCpBZXkdw/HDSIhEhuuJ+QjdEjRNpKpOiGBvaxqFMNMc40VIi0mUhyAaLAOhZp\nx2inGhAJkVxxEGnQKEaTEAmR3LCN5atDqhBpAJEQqXIRKdoFqw2IhEhuzI9lmLGL0SNEQiRH\nZscyTNhF2SEhEiI54iDScGS3RnH2Ji2RehDJA5axSMvsorzn90BaIvl0A5Hm7SViPouESIjk\niotIka5qaEAkRHLDKhZ5p5j7I0RCJFdmihTxMrsWREIkN+bFEumdGj4gEiK5MVek7jPRmoRI\niOTGrFiGrijeLgmREMkNi1i+1n1XiKSASIhUOYgU67rvDkRCJDfMsfxciMRkgwIiIVI1X6So\n5+wQCZFcmSFSpA/pU0AkRHJjRizRz31XiIRIrtjHMtwQMrUvH5QAAA7QSURBVGaPEAmRHJkh\nUiVivoCiA5EQyQ1TLCMXUKxboH1BJERywxCLtDn2me8WREIkN6xFSmDuu0IkRHJlOhahvozz\nCRQyiIRIbtiKNPRE9EjfIBIiVfbfESLpQCREqmaI1O8ZU+wjIBIiuWEdSwInkSpEQiRXJmOR\nN6YwZ4dIiDTG7SREfp/eZyqWr23xa4RIiKTQtfhzdwK1mN7VcVukIBIiSbQiFaIoq+pViNvk\nrk6bogWREEmiFSkTZfO6FKfJXZ02RUtUIgkFZQsiWaFMsE2PbIKPZVviEklpuvotiKSjdecy\niJRN7rpBcQ4EIiGShBD59XYXf/XLspiebQg+lm1BJESSkA6LhcjKyV1tN6Qw+Y1IiKTyfN5u\ned5OORSTHmm/o59zSJVptBUFiIRIbtiJJIYnMMcOIiGSG5pY1Gp/X9MXU+SjxC6S6bGxG4sU\n/DXXZdFM1V1PQpz/pve0Fqld/x1yzF6IXqTV3HASSSpZkLyyuuWXWef7eXJXiyjed/xOYJCE\nSIgkcRF5Wf9zedVOXRZPfwvpf7GDSIgkIZrVQWJYIrTwhOznERTxd0iIhEgy/Vo76Rf9rqZ3\nh2O6JDxCJESSuYhnVV2bf5oeaXKQNBrF2AV9iKQBkeIV6Smy4lnlWW3S/SQmL+0zidSdPQp9\nltIbiIRIMvfss3z++rNVu7p+2K6+TmVVQwMiIZLK3+XUaJJfX9P7WYhUJXJY14BIiOSGMQph\nt1skIBIiuWEWKaEDO0RCJFfsVjasX45AQCRE0jH3PNJR4loFREIkHYg0A0RCJDd+oogjLFcQ\naQeRwr+awgJEUkCkHUSSC3hYEEkBkRBJ4XHN2+4yLx7TO5qKf/gudx6IhEgS5UlaA7Tkwr60\nTiJViIRICoXI/tql39Xrni25sC+tZQ0VIiGSQtZdQdHynHVhnxj9LcQY1wGREElCORibcx5J\no1WIMa7D0UVS75uvtFPprTBFCnAS3LlHQqT5HwlLpDVa+EYidf+uWI3zqcdI9+7yiXljJM1c\neFjBrQoiIZLMWerfT/b3/v4VSQjTsWFcIBIiKTyK9jxSll8XnEfqLAottFVBJERyY1Ik4x7R\ngUiI5MZEwdObakAkRHJFaF5XiGQJIiFSpXxH3zEgkhWIhEjVpEiMkaxAJESq5O9o5FrZ1Jas\nIhIiuTIhUoBLNlYHkRDJDW3BE+yOKkRCJFc0BR+6osPG5QgiIZIbowWXHs932MDcQCREcmN0\njlsM/yGSmb1FUm7Cg0h7MSaS+Pp/QhxRpH9SC0SkvRDKD+kXqVdKCERCJDf0IqXxGPMvEAmR\n3BgTaVjSkJxGiIRIruhm7RI8h9SASIjkhuE8UmogEiK5cdiCrwMiIZIbKZ4smgCREMmNBM8V\nTYFIiOQGIikgEiK5keBJ1ykQCZHcQCQFREIkNw5b8HVAJERy47AFXwdEQiQ3DlvwdTiGSCPP\nnNA/beKQIoX4aIppxNi3tEtJguAgIu3RwrcVSf7lEHxbk+4yuxZEQiQ3fkSS/k2QTUVyfmI2\nIoXHu6zdEd3YBbMpsaFIC56YjUjh8R7cVf21fMrbybGhSAuemI1I4aFc2IdIG4rk/sRsRAoQ\ndUwkGCPNxlUkwxOzlRnu8XcTw7GeN+JHJGbtZrJBjwTh8y0S55HmsmCMZPvEbAifrzFS6mw5\n/W3/xOyAoIlo+Jq1S51tzyPZPjE7IGgiGr7OI6VOmCsbAoJGooGKUUAkA7QXDVSMAiIZoL1o\noGIUEAncQCQFRAI3EEkBkcANRFJAJAO0Fw1UjAIiGaC9aKBiFBDJAO1FAxWjgEgGaC8aqBgF\nRDJAe9FAxSggEriBSAqIBG4gkgIigRuIpIBIBmgvGqgYBUQyQHvRQMUoIJIB2osGKkYBkQzQ\nXjRQMQohi7TLPbCCYquadmHvugmN+TW4nUi+dvKWUIC5BYJtUa1D2i1B/znrQKSQcgsERJoP\nIoWUWyAg0nwQKaTcAgGR5oNIIeUWCIg0H0QKKbdAQKT5IFJIuQUCIs0HkULKLRAQaT6IFFJu\ngYBI80GkkHILBESaDyKFlFsgINJ8ECmk3AIBkeaT1upvgJVAJAAPIBKABxAJwAOIBOABRALw\nACIBeACRADyASAAeQCQADyASgAcQCcADiATgAUQC8AAiAXgAkQA8sJZIRSayopTeeF6EuLzG\nt81ISLnN+ayEfvYupTeWFKl943zXbDsiZfNdPS12vJ2so70ZL52bUXXmxLrd7IpnHe4UK4l0\nbhv76fPGvX0jK8e2zUlo8Cibm9DP3q+sS+i1tEj9G9fRbYekqxlz0yo+X6qJp/EJDzOqzpxY\ni23xbMOdZB2RHiJ7Vs9MPN7vZPUbZS6KsW2zEmq5N2/MSuh370tTmLqyL0uLdBPnsvmr9pxb\npFBp66QQuWm/p7iUTfgXc5J1nRja/oyqMyc2q3i24U6zjkiFaI50/rq/0lX7smm1ZdOR/Gyb\nlVBLmeXabdYp9V9G82NZkc7tl/9qIpxXpFDJRPNH3NxY808NGqj/1pj2sq86i8RmFc823GnW\nESkXzQHT82P55dNz/myblVD/bjk3od+9s76es6VFGow8zy1S2LRHz1Y7mptg/UfGtJd91Vkk\npuxut691uBrWEUl8/zE4ieqatV3t77ZZCTU82+5tXkK/e1/7Q7vr0iJJb8wrUtAU4ma3Y9n8\nBTHwNNeJfdVZJCZhU7xqRrg6NhJJiHyYIlguUtchLRWpujWDzOy2uEin9m/pIyqR/kT3x8qC\nm7jb7OZPJOu9WqyKNyNcHZuJ1Ew2XJb/+a/aQaR225yUru/JtmVFuoq8rJ7nqES65ZnlQO+V\n2R3H7iSSXfHsw9WymUjNGOnVzG4uFqno/8YsFOnW/BGq3b4tLlI7fZpHJVLVDGttDnbKzOrI\naS+RrItnGa6edUTKfkV6//jZNiuhz3vzEvrd+9QeIJaN2wuLVNuYXefHFhzqA71L7fBb3u88\ndeJH3tFUJ7OqzrqCJ4unoA/XjjVn7V6fKRhpKvJn26yEpImdWQn97i38Fakr1mlukUJDFUnf\nWj/7vU7nl2WCdrN2llVnKZKheE5p6lhHpGt79HX/jOC6N17NDMrPtlkJNYdkN+22OSl1fwLb\nP0TLitSdh7g1TWBekUKlC+hlXmVwt5sRazG101lVZ9foLYtnHe4kG61sqIvZnv3/W76yIR9O\nSS1c2VCIZhlWsXyxRbc44jQ7tmBpAypz46DhNcMjY9ufVXVWItkWzzbcaVZaa3dqO/U2kC7o\n6+cNaZtDQv3YZnZCvymdPRWp7NZq5fOLFCqZXRQXeQGxCeNec6rOKkvr4lmGO81KInUrq9uX\nfSD38/CGtM0lIWk4PCehkZQ+bywr0qv+yvK7Q5GCpY7iZP4DLbyKNKfqrLK0L55duNNwPRKA\nBxAJwAOIBOABRALwACIBeACRADyASAAeQCQADyASgAcQCcADiATgAUQC8AAiAXgAkQA8gEgA\nHkAkAA8gEoAHEAnAA4gE4AFEAvAAIgF4AJEAPIBIAB5AJAAPIBKABxAJwAOIBOABRALwACIB\neACRADyASAAeQCQADyASgAeSEKl7ml617LHVaVAWJyFORWnes+rrtXki3owngpcX8Xngcpvb\n+fZJbIrJTD4bhbh/v2XJoueapyDSqamgV1mLVNo/LT5N/obHRdo8CbKt17ki5XXi1/51//Rd\nkZVDYlNYi5RZ7D87AxMpiNRW0Fnkp/zwj0lemdqjov5b8yqsTJIbr3UWQnz+mF3Eucnt3HRR\n5iSsRepFRSTfdBX0vIjLw+6QJVXqHqI/LLoLYa4qR5Hk120mpV2nZi3SqXMVkTzTP9n6nt2E\n8Tg8bW6f0UvR/FkX72O36l4fknXPHK+7lFxk13e9fiy4nUTWj3fO9dBHrux6U/vccPUx41LL\nHd4fyadqHzte/JajPIlc2din9Gzf/hSqe2J5v3/97rVNtujHanKKC+ouHZGul+pyNe+dMrl4\nDi8fzVHwR6RrZ0DR/tqObK4/IuXtG83B8+17nHUeNqkiFeIyHOcNX9JIPv3H85/tefvis3FI\nqT5kfFTvI/qhUP3+dYrNO/dzn46S4oK6S0CkZRWUEuphlyySEH/tCKr99VzWqpzUzc3RYP1+\neW56/awR8q/ZpeNPZM/qmTVpqF/GuZkhfEh5j+YzfHykHMrGd8HLNuvmLSnnfv8+2e7f7CvF\nBXWHSPBGL5K0Xbz/3Ksi5e2Ap2wOq8TXMXTe/n7vOwZ5y/3STNrdR/KW8snbV/fRcoxvvDW9\nYffRd879/kOyr7HI5lWXAiLBmymRXvfrWWpuvyINE+eiHX7kz+dPuuOt9XHNmrY9vK/NZ7wc\nXz+GV6fa6q+PKjtOROYIIsEbaYz07HqW5qU01rAUqbo2w5vsPc89KVKT1/s4cSKf8XJoRHqI\ni61IPym6gUjwpp+1e76aTuUuN7eLON3uL4NIclL34vQZI2lEer/8fHwin/FyaERq/ibYifSb\nohuIBG/680i5yP/ew/X6j/vQCqdFyr9PLnyqfRip5Mq79dvdvF7ZDfrfn/nOp/v4eDnkjXK2\nL3H6fPT+6V5/RPpN0Y00RGJlkB33bmXDtT8uO9UtvTx3ze1RPX/HLq/q0xrbCbK6U8ubz/3Z\nzNrVrf9W1j/O3dRAl9hIPnd51u5ru7yxo391Fe9CDbN20mY5iO8U3UhBpNN7+RUYuL9HOk1n\ncXufoin6dx9yG+zq9f1nvRtsNAb+vXfu+ZzNUVrrkGqzoUtsNJ/uDNVltBzyxo7hVSYVSpou\n/BZpJEUnUhDpcUIkW/rV3/dzexx2zfpBez2UEOfHXZ2B6Or1c5B3qz/anWFtVzY8pGRvWb++\nQG2tz0vd3M9NRzV8SWP5tOUoxsuhbFQy6GfE3znrJhtGUnQhBZHAgTurQGaBSAAeQCQADyAS\ngAcQCcADiATgAUQC8AAiAXgAkQA8gEgAHkAkAA8gEoAHEAnAA4gE4AFEAvAAIgF4AJEAPIBI\nAB5AJAAPIBKABxAJwAOIBOABRALwACIBeACRADyASAAeQCQADyASgAcQCcADiATgAUQC8MD/\nKJy2DFry7toAAAAASUVORK5CYII=",
      "text/plain": [
       "Plot with title \"Histogram of t\""
      ]
     },
     "metadata": {
      "image/png": {
       "height": 420,
       "width": 420
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "86536f34-ff89-4b67-b4e7-fc2496f101f3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS\n",
       "Based on 1000 bootstrap replicates\n",
       "\n",
       "CALL : \n",
       "boot.ci(boot.out = results, type = c(\"perc\", \"bca\"))\n",
       "\n",
       "Intervals : \n",
       "Level     Percentile            BCa          \n",
       "95%   ( 0.6753,  0.8835 )   ( 0.6344,  0.8561 )  \n",
       "Calculations and Intervals on Original Scale\n",
       "Some BCa intervals may be unstable"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "boot.ci(results,type=c(\"perc\",\"bca\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d6e4153-b856-49f9-b571-9e34d41fdae8",
   "metadata": {},
   "source": [
    "- Level：表示置信水平，这里是 95%。意味着如果我们重复进行抽样和计算置信区间的过程，大约有 95% 的置信区间会包含真实的总体参数。\n",
    "- Percentile：百分位数法计算得到的置信区间为 (0.6753, 0.8835)。该方法直接使用自助法重抽样得到的统计量的百分位数来确定置信区间。例如，对于 95% 的置信区间，通常取第 2.5 百分位数和第 97.5 百分位数作为区间的下限和上限。\n",
    "- BCa：偏差校正加速法计算得到的置信区间为 (0.6344, 0.8561)。BCa 法是一种更为精确的自助法置信区间计算方法，它考虑了自助法统计量的偏差和偏度，对百分位数法进行了校正，通常能得到更准确的置信区间估计。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "82bc3d06-17bc-4811-856b-7e2a5175a8b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "bs <- function(formula,data,indices){\n",
    "    d <- data[indices,]\n",
    "    fit <- lm(formula,data=d)\n",
    "    return (coef(fit))\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "13c0a738-989e-45a5-a69f-96885ad0bd51",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "ORDINARY NONPARAMETRIC BOOTSTRAP\n",
      "\n",
      "\n",
      "Call:\n",
      "boot(data = mtcars, statistic = bs, R = 1000, formula = mpg ~ \n",
      "    wt + disp)\n",
      "\n",
      "\n",
      "Bootstrap Statistics :\n",
      "       original        bias    std. error\n",
      "t1* 34.96055404  0.1627187591 2.471989763\n",
      "t2* -3.35082533 -0.0815507218 1.144253931\n",
      "t3* -0.01772474  0.0002047975 0.008758798\n"
     ]
    }
   ],
   "source": [
    "results <- boot(data=mtcars,statistic=bs,R=1000,formula=mpg~wt+disp)\n",
    "print(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "801090fe-89f6-40c1-9383-329332a5690e",
   "metadata": {},
   "source": [
    "三个回归系数：截距项、车重和发动机排量。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "23057a26-fdf9-450d-8836-372ce963e42d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS\n",
       "Based on 1000 bootstrap replicates\n",
       "\n",
       "CALL : \n",
       "boot.ci(boot.out = results, type = \"bca\")\n",
       "\n",
       "Intervals : \n",
       "Level       BCa          \n",
       "95%   (28.90, 39.26 )  \n",
       "Calculations and Intervals on Original Scale"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "boot.ci(results,type=\"bca\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "5bf12436-e1d1-4ab6-9199-c84608ba6e33",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS\n",
       "Based on 1000 bootstrap replicates\n",
       "\n",
       "CALL : \n",
       "boot.ci(boot.out = results, type = \"bca\", index = 1)\n",
       "\n",
       "Intervals : \n",
       "Level       BCa          \n",
       "95%   (28.90, 39.26 )  \n",
       "Calculations and Intervals on Original Scale"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "boot.ci(results,type=\"bca\",index=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "7612bc0d-ef29-4490-8aa4-3cd6e836d608",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS\n",
       "Based on 1000 bootstrap replicates\n",
       "\n",
       "CALL : \n",
       "boot.ci(boot.out = results, type = \"bca\", index = 3)\n",
       "\n",
       "Intervals : \n",
       "Level       BCa          \n",
       "95%   (-0.0343, -0.0006 )  \n",
       "Calculations and Intervals on Original Scale"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "boot.ci(boot.out=results,type=\"bca\",index=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "3d242374-e39d-4114-8a87-98629c914738",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAANlBMVEX9/v0AAABMTUxnaGd7\ne3uLjIuZmpmmpqaxsrG7vLvFxsXOz87T09PX2Nff4N/n6Ofu7+79/v03sFGxAAAAEnRSTlP/\n/////////////////////wDiv78SAAAACXBIWXMAABJ0AAASdAHeZh94AAAgAElEQVR4nO2d\niZajqhZAH2aqdG4l8f9/9sUxUlECeESGvVevrgzK4RB3UETzvxoAVvO/vSsAkAOIBCAAIgEI\ngEgAAiASgACIBCAAIgEIgEgAAiASgACIBCAAIgEIgEgAAiASgACIBCAAIgEIgEgAAiASgACI\nBCAAIgEIgEgAAiASgACIBCAAIgEIgEgAAiASgACIBCAAIgEIgEgAAiASgACIBCAAIgEIgEgA\nAiASgACIBCAAIgEIgEgAAiASgACIBCAAIgEIgEgAAiASgADxiqSU0h+9X5hyDlKZn0qpaaQw\nUSEdEhfpt5qTS5yfV+yJSIGiQkIkLtJ8LyXOQal7+KiQEAmJZF4oUF1CRoWESEik/s/z5/h6\ndPrXvaKGxW7nZu/r1q/yeD07XidrPg7q8nr07/R6fLg8hvKuB3X4retrpY6/enitPM2caVSA\njuREelT9hnzUNulj//jUrvHbL/Je89CuMCylfrtXu+ePy/jayLS8P+YgEnySnEivjuLVGT1f\nG/p1skmfBkM6k6rx6bCmala7vtR51vVF80Opaupgh1YeIsFXohbpQ4fu/2bH7PnqZMbX6tvr\n7/X52ut7/X3tjf17ydH8qd5rNgI1YwYPraTXq9emu7q3f96x/5THMRJ8IzmRGjnGQ6Fhkz43\n/VPDpR2lPnWbfyPUsObtT9Hd/7/an/cCf8pDJPhGciL9dC/0Lr3ferbPH+0L1bCh/337tcC/\ny1GNItUff8b1puUhEnwjapH0R/2fy3Bk8/h4a3ikPkXqnv87TMw0i6Q9QiT4Qnoi1c9/3ZDa\nsZ7tkarZHql92uzqHc7Xu1OPVP19E5HgkwRFamjP8rxfO309RmrfPfSvfxXpxDESuJGcSIe+\ns3h3Fc/FUTv1R5L+7/ceyWLU7lkDvElOpNc2fny0Yw7NTIVmDK/5O55p7aaWfp5Hags6tgvf\nqq8i/S3vj0hDVICB5EQaBxuaQ6RmnLp7cJxu922PorSZDe3Lv8M4RTvibRLpT3l/3hyjAvSk\nJ1J3fHTsjmGao5nenXM1OcF0b+ba3T4kaV6uzvfHMGFhpvQBrby/b570q5MA4hVJgGd3IAWw\nOVmKpNr5ePX9qE+gA9iMLEV6DxXoU4MAtiJLkcZLLRhbg0BkKVL9/Gmug6jO9EcQiDxFAggM\nIgEIgEgAAiASgACIBCAAIgEIgEgAAiASgACIBCAAIgEIgEgAAiASgACIBCAAIgEIgEgAAiAS\ngACIBCAAIgEIgEgAAiASgACIBCAAIgEIgEgAAiASgACIBCAAIgEIgEgAAiASgACIBCAAIgEI\ngEgAAiASgACIBCAAIgEIUKJIau8KxIcCDfcWRCSoaZI/IJIVbDUf0CQaiGQFW80HNIkGIoEf\niNSi/vx1AJGgRqQeRIJ1IFKD+nhgT4kisdV8QJM0IJIbbDUf0CQNiOQGW80HNIkGIlnBVvMB\nTaKBSFaw1XxAk2ggEviBSFobIBL4gUhaEyAS+IFIiOQOW80HNIlafGIHIkFNkyCSD2w1H9Ak\nGohkBVvNBzSJBiJZwVbzAU3SXm4/PnRfu0SR4IPSRVJ1a9GgEiKBH4WLpIYWQCRYRfEi6Zf0\nIdISa2+2lDtlN4mqEckS9V9L96fsrWaWspsEkaxBJDM0CcdIViCSGZqEUTsrEMkMTcJ5JCsQ\nyUzJTTKTOyItgUhmCm6SudQRaQlEMlNwkyCSC4hkptwmmc0ckZZAJDPlNgkiOYFIZmgSDURa\nApHM0CQaiLQEIpmhSTQQaYk5kZjB+qbUVljIG5GWmBWJ7mmk0FZYShuRlkAkM4W2AiK5gkhm\nymyFxawRaQlEMlNmK7RZzx0mI9ISiGSmzFYYBps+VEKkJRDJTImtMHEIkWxBJDMltkLr0fy1\n5Yi0BCKZKa8V+p264UDpz5vu5SHSvhWLhNJaQan2AAmRXEEkM6W1gnr/6zumv2+7gkj7ViwS\nCmuFYZjhz52KtfcdQaR9KxYJhbVCZ4+qVT0/3RKRlkAkM4W1wnvasvF9FxBp34pFQmGtoOZP\nH73fdy8SkfatWCQU1QrjkB0iOYNIZopqhf4IyZAzIi2BSGZKaoVxhp1hEfdSEWnfikVCQa0w\nDngjkgeIZKagVlDDaVhE8gCRzJTTCqNGppQRaQlEMlNOK6h6YTbDn4VcQaR9KxYJ5bTCcFWs\nMWNEWgKRzBTUCrMXIM0u4wQi7VuxSCioFdrO6MvtDBFpCUQyU1IrWNwVFJGWQCQz5bSC1b11\nEWkJg0jcuLguRiT1ZY7duJx70YhEv1QXIlJvkfqeLSItgUhmishf2cwO6pZ0LxyREKkuQ6Sh\nK7LokhBpCUQyU0L+iCQAIpkpIf/h/JHF0BIiLYFIZorIvxuwsxmjRaQldJF6EGmkiPzHUbvv\nS7oXXqRIn0+K2JAMFJK/7RlDRFoCkcwUkr9tmoi0BCKZKT3/PyDSEohkpvT8/4BISyCSmSLy\nt59SiUhLIJKZAvK3nK/aLetePCIhUl2CSMryXGy3sHv5iIRIdf4ived9M/y9CkQyk3v+3Qmk\n79eYD4u7R0AkRKqzF2mY0GB7kIRISyCSmczzd9yzQ6RFEMlM5vkrl4l2NSItg0hm8s6/+41L\n65l2iLQMIpnJO39Vu3RHNSItg0hm8s5/uFMUIq0Gkcxknf9w/Zl9koi0BCKZyTl/hUhyIJKZ\njPMf5gYhkgSIZCbj/N+nYhFpPYhkJt/8h56IHkkERDKTb/7DtRNO93dHpCUQyUy++Q/nYRFJ\nBEQyk2/+w9QgpwzDi3Q9KHW6rSsjBIhkJt/8nYe+25Xc4/iK1NXs2FXy4llIOAoV6ffn1H5A\np8uvecFM8x/27BwmNfRrObNKpIu6POv6cVFXz1KCUaRIz8P7XKQ6GhfNMv8GNUxYdVvLPdAq\nkSr1bB4/1cGzlGAUKdJFVf/u7aPHrTLvNmSZf8NwaWzcIg17nvH/cmSRIlXqPj6+q8q0aJb5\nTy5CcswvtEjnQSTjhxQDRYqkfb+Zv+yyzL8efzPW9Zs+qEinn+tN/Xs9fF7iH20oUqTie6TB\norhFGiuoVPX0LCUYRYr0Oka6PdpHhR4jDQN2zoceIc8j3e/X6+nUDjlcoveoTJGGsxMtB+OH\nlGX+k+xd13QPxsyGjEWqfy/teaTq9FPieaTxAMk5O0RaolCRrMkx//cBEiKJgUifKOW/75MC\n/h0SIi1SuEhft6UM8x++IRBJEkT6skCQagRlRYe0x/B3InsGRYrk8Bnll/+KI6SgIl0RKXp+\nq5JFGk/Heq3sjP95pMo8nzguihSpfp7UsT0jW96u3Xh85JVZ0GOke/wTg96UKVJd/1PtLK7i\nRFq5sxR2sOE6mcr1UZPI9vtKFal+HNXpWaJIbnfN/7u6+yoBRu1i+JCKFamuf1R1K02ktd/i\niLREwSLV98P37Smz/BFpK0oWqblurFiRPNd3XwWRChDpO5nln6pI8Z+jQCQzmeW/auy7RqRl\nEMlMXvn7Xhj7LsB9FXbtEKmO4zOSYzAIkcRBJDN55T/cNcg7K0RaApHM5JV/v1PnPxUgqEhJ\n3Q4Xkcxklv/aKTUBRUrsdriIZCav/FfPTQsoUmK3w0UkM3nlv3qOZ0CRErv5ICKZySr/8RIK\n/xLcV1n3sy5zTz4X9QwhCSKZySr/4acukxCJHikrssq/c8j5zvnTEtxXWXGMlNLtcBHJTE75\ndxp53Kh4UoT7Kt7D32ndDheRzOSUf98hJbJrl9jtcBHJTEb5S1yZzcyGJRDJTEb5DxYh0hYg\nkpl88he5VQgiLYFIZvLJf+2lSF0h7qsgEiLVcXxGIqy+FKkrxX0VREKkOo7PSASRPTtEWgSR\nzOSSvxLZs0OkRRDJTC75y3RIiLQIIpnJJX+RIyREWgaRzGSSv1CHhEiLIJKZTPIfj5AQaSMQ\nyUwe+TcWdT+LtLYg91UQCZHqOD6j9Qx7dat/AAWRlkAkM1nkr6T27BBpEUQyk0X+UkMNiLQM\nIpnJI//1N2sYCnJfBZEQqY7jM1qN2J4dIi2CSGayyH/csUOkzUAkMxnkP7n1ASJtBiKZST9/\nVSNSABDJTPr5S3ZIiDQXvQORjCSf/1siRNoo+pI7iDQh9fyVaIeESHPREcmC1POX9QiR5qIj\nkgWJ5y/cISHSXHREsiDx/IU9QqS56IhkQdr5S3uESHPREcmCpPOX3rFDpNnoiGRB0vmLe4RI\nc9ERyYKU8xfXCJFmoyOSBSnnL+8RIs1FRyQLEs5/gw4JkeaiI5IFCee/gUeINBcdkSxIOH9E\nCgMi2ZBu/lt4hEhz0RHJgnTz38IjRJqLjkgWJJv/Jh0SIs1FRyQLks2/VQiRAoBINqSa/xZj\n3zUizUZHJAsSzX8jjxBpLjoiWZBo/ogUDkuRtvg4EiLNxLfyCJHmotuJ9F/ZHVOSicvP+h5L\ndl8FkRCp3vsz8mQzjxBpLjoiWZBg4tv1R4g0Gx2RLEgvcYVIYUEkG9JLfEuPEGkuOiJZkFzi\nm3ZIiDQXHZEsSC7xTT1CpLnoiGRBaolv6xEizUVHJAtSSxyRgoNINqSW+KYaIdJsdESyILHE\nN/YIkeaiI5IFaSW+7X5djUiz0RHJgrQS39ojRJqLjkgWJJX45h0SIs1FRyQLkkp8c48QaS46\nIlmQUuLbd0iINBcdkSxIKfFOoS09QqS56IhkQUKJB+iQEGkuOiJZkFDiATxCpLnoiGRBQokj\n0j4gkg3pJB7CI0Sai45IFqSTOCLtBCLZkEziIYYaEGk2OiJZkEziQTxCpLnoiGRBKomH6ZAQ\naS46IlmQSuJhPEKkueiIZEEqiY+TGhApNIhkQyKJB+qQEGkuOiJZkEjiiLQfiGRDGomH8giR\n5qIjkgVpJI5IO4JINiSReIDrJ4ZI7qsgEiLVe39GlgTrkBBpLjoiWZBE4oi0J4hkQwqJh/MI\nkeaiI5IFCSQeaHZQF8t9FURCpHrvz8iKgB4h0lx0RLIg/sRDdkiINBcdkSyIP3FEQqQEiD7x\noB4h0lx0RLIg9sQVIu39ISGSDbEnHtYjRJqLjkgWRJ544A4JkeaiI5IFkSce2CNEmouOSBbE\nnXjoDgmR5qIjkgVxJx7aI0Sai+4kUshPKybiTjm0R4g0F91JpFL7pahTDu4RIs1FRyQLYk45\n+I4dIs1GRyQLYk4ZkXoQKX4iTnkHjxBpLjoiWRBxyog0gEjxE2/Ke3iESHPREcmCeFPuFAp8\nUgKRZqIjkgXxprxHh4RIc9ERyYJoU97FI0Sai45IFkSbMiK9QaT4iTXlXYYaEGk2OiJZEGvK\n+3iESHPREcmCWFNGpAmIFD+RpryTR4g0Fx2RLIgz5Z2OkBBpNjoiWRBnynt5hEhz0RHJgihT\n3q1DQqS56IhkQZQpI5IOIsVPjCnv5xEizUVHJAtiTBmR/oBI8RNhyjt6hEhz0RHJgghT3tEj\nRJqLjkgWRJgyIv0FkeInvpT39AiR5qIjkgXRpbznERIizUZHJAuiS3lXjxBpLjoiWRBbyvt2\nSIg0Fx2RLIgt5X09QqS56IhkQWwpI9IMiBQ/kaW8s0eINBcdkSyILGVEmgOR4ieulHceakCk\n2eiIZEFcKe/tESLNRUckC6JKefcOCZHmoiOSBVGlvLtHiDQXHZEsiCnl/TuksCL9/pzaZE+X\nX/OCiBQ/MaW8v0chRXoeJl8cR+laCYJINkSUcgQdUkiRLqr6d28fPW6VupgWRaT4iSflGDwK\nKVKl7uPju6pMiyJS/ESTsipNJC1Nc857nVTrQSQLokk5Co/okbSwX9xBpAmxpBxHhxT4GOn2\naB/FeoyESA7EknIcHgUd/j5Ocj48hWslASI5EEvKBYpU/17a80jV6SfO80iI5EAkKUfiETMb\ntLCIZE8kKUfiESJpYRHJnkhSjsQjRNLCIpI9caQcS4eESFpYRLInjpQRyQgi7cv1oNTpZl4m\nipRjGWoIPLNBQ7hWEiBS97H05ymMp/piE2n3mriv4ivSFZGip/1YLuryrOvHRV2Niwaqkol4\nOqSgu3b3ynzxxBtE2ol2g6xUe7b8qQ7GRYNUyEhEHoU9Rrp/2VsYSVWkWD5Vb9q6KzV5srxo\ngOp8oViRXnt398X3rPf7tmO9SKl3T23Ln4dsYpxYPK3B3tuLVhn3VRi1y1mk08/1pv69Hj4v\nUU4sntYAkb6BSDsx2TCVqmKcWDypQEwerRXp8PMQq8oERNqL+/16PZ3aIYeL0SNE0lkp0iuJ\nLVxCpPjZO8u4PFor0vPf2c+lKEeEEOkLEQwIzVRl34r0CBwj/f4c3F1CpKi5VupgPB27d48U\nj9E9MoMN9+qVz5eGdwGR9uJ+UtW1/mk30KjvPRiZRzIi3Y4WDe8CIu3Evf0gL+r8rB+niKcI\nRdchCYj0/Hl1R4fb82XTSaZOiLQb5+bc0aU7ExvzFKHoPFot0m8z2HDp5it8TSr2e38jUvcR\n9t+IUR7HdrGzE6kZZrgO5xvMU0oSuPc3InXb5b9uny7eKULxebT6PNK3y78mxH/vb0Q6N0dH\nHc9zlJ9RFzs6j1afR3JYr4A7rSYv0rMat031ZQdjxywj9Gj9zIb+QWVudW3ZjycStZIAkV77\nDcMHWX255AWRNIREeljkRI+UFYiksUKk23TsxDxU2lLAvb8RKUjkCD1a1SNNR+EOXwa0G/K/\n9zciBYmcm0i1xakjjezv/Y1IIQLnKNJGIFL87JVlnB6tEalJZKOsECl+EEkDkaZhEcmevT6j\nOD1i104Li0j2IJIGIk3DIpI9+2QZq0erRboe6vpxsBr9dgCR4geRNFaKdGvyaS6PVaImIVL8\n7JJltB6tFemo/tV3daj/CV4eW6coUg8ibRw0Vo8kTsi2d/QufdRupnvKnD2yjNcjCZFO6oZI\niBQiZL4iHdX91kzkLn7XDpFChMxXpHYG+E+Tov2VshYgkhSPo8u1ly6EzzJmj9YPf3cXRBz+\nCdWnA5GkeKrf288mJSOSBidkp2HzE6k+ndV5k4KDZxm1R4ikhc1OJKUOknfA1YreptjleFF7\nhEha2OxEanqkKoseSWUu0s9hi+wQSYqnuudxjBS5R2tF+tkmPUSS4iF6WmJK2Cxj75DWilRt\nsweOSPGzm0hB41oT9J4N1iCSIC73HnQpVrKwr8Fi92itSCe1yfk+RBLE5d6DLsVKFvY1WO4i\nPaqj7JVIHXudo8hOJNd7Dzqwk0gho7qwetcuj8EGO3fSE8n53oP2hMwyfo8QqQ+Yq0h1Fsex\n+Yu0EYgUPwGzTMAjROoDIpIriKSxWqTbqb247yFUnw5Eih9E0lgr0rFLT1WiJiFS/ITLMgWP\n1op0Vcdnk99Vdq4+IsUPImmsniL07EaFGLVDpK0ClSBSd//vGpEQabtAKXi0VqRD3yPdEz5r\n3gZEJFdCZZlGhyR0jHQTngWOSPETKMtEPFo9anfqk5S97AWR4geRNETOI6mT7E2EECkBwmSp\nihFpExApfoJkmYxHiNQHRCRXEEljnUi3czNN/3iRviYJkeInRJbpeLRKpMdxzPLIXDtE2iBG\nMh6tEelZqcOtudL88e+gRG8HgEgJECDLhDxaI9JlMuZ9VKI3T0Ok+Akq0vax1rJCpIN67889\nEv9ZF0RyBpE0VoikpcdcO0QSj5CQR4jUB0QkVxBJA5G6gIjkSkCRNo8kACJ1ARHJlc2zTKpD\nWieSxs61WhkQkVxBJA1E6gIikitbZ5nSSaSauXZDQERyJZxIGweSAZG6gIjkysZZJtYhIVIf\nEJFcCSbStnGkQKQuICK5sm2WqXVIiNQHRCRXQom0aRg5EKkLiEiubJplch0SIvUBEcmVQCJt\nGUUSROoCIpIrW2aZXoeESH1ARHIljEgbBpEFkbqAiOTKhlkm2CEhUh8QkVwJItJ2MaRBpC4g\nIrmCSBqI1AVEJFdCiLRdCHEQqQsoKlJ6m4EH26WXYoeESH1AUZH+K6FjQiQNROoCIpIriKSB\nSF1ARHJls/SS9AiR+oCI5AoiaSBSFxCRXNkqvRTPxtaINAREJFe2F2mjANuASF1ARHJlo/QS\n7ZAQqQ+ISK5sLtI25W8FInUBEcmVbdJLtUNCpD4gIrmytUibFL8diNQFRCRXEEkDkbqAiOTK\nJukl6xEi9QERyRVE0kCkLiAiubJFeskONSDSEBCRXNlWpA0K3xZE6gIikiuIpIFIXUBEcmWD\n9BL2CJH6gIjkCiJpIFIXEJFc2VIk+aI3B5G6gIjkinx6KXdIiNQHRCRXEEkDkbqAiOTKZiLV\niCQGIsWPeHpJd0iI1AdEJFcQSQORuoCI5Ip0eu0+XbIeIVIfEJFckRcp5Q4JkfqAiOQKImkg\nUhcQkVzZTCThcgOBSF1ARHJFOL3EOyRE6gMikiuIpIFIXUBEckU2vf5MLCLJgkjxIy+Sav+I\nFhsOROoCIpIrG4iUcIeESP3Hh0iubCSSaKkBQSQnNxBpRDS95DskREIkTyTTS98jREIkT2RF\nSnqeXQMiIZIfgukpeqSNQKT4QSQNREIkPzYRSa7M0CASIvmBSBqIhEh+yKWXg0eIhEieIJIG\nIiGSH4ikgUiI5McGIomVuAOIhEh+IJIGIiGSH2LpZbFnh0iI5ImcSMlPD2pAJETyA5E0EAmR\n/JDftZMqcBcQCZH8kBQp5WvMexAJkfygR9JAJETyQyq9/lqkNH/M5Q0iIZIfciKpwaWUQSRE\n8kNUpOQ9QiRE8kRSpPQ7JERCJE+E0stjXgMiIZIvsiKlPtaASIjkiUx6Y3eUemshEiL5IZKe\nymXPDpEQyRM5kTIYakAkRPJFIr18OiREQiRPpEQauqXECS/S9aDU6WZeBpHiR65HSn/Gah1U\npK61jl3jXcyLeobwAJE8kd21EyhsX0KLdFGXZ10/LuoqXCtfEMkTUZEEytqZ0CJV6tk8fqqD\ncVHPEB4gkifr08toqCG8SEObmdsOkeJndXoKkTxDNbHOg0iVcVHPEB4gkieIpBFUpNPP9ab+\nvR4+L+bRBkSKn7XpvUfscvAorEjjt49S1VO4Vr4gkicyIuXSIQU9j3S/X6+nUzvkcDF6hEgJ\nILhrJ1GdvWFmAyL5IdgjidRnZxAJkfyQE0mkOnsTVKTfn1PbcqfLr3lBRIqflell5lFIkZ6H\nSW9+lK6VL4jkybr0MtuxCyrSRVX/7u2jx61i+Dt1BETKZuy7DipSpe7j4zsnZFNnVXqjRbl4\nFH7299yTz0U9Q3iASJ6sFWlwSao+O0OPhEh+CB0jCdVmd8IeI90e7SOOkTIAkTRCDn8fJ0M1\nh4+pDdMpI4gUPwK7dojkye+lPY9UnX44j5Q8K9J7W4RI24JI8eOf3vt3XPLxCJEQyZMVImV3\nNrZGJETyxTu996nYjDzaTSTOI6XOGpHq7DokREIkT1b1SHmdjG1g1w6RNLafoT+MM2TlESIh\n0pQQM/Rz3LFDJETS2H6GvspvUkMLF/Yh0oTN50OqOsv+iAv7EElj8xn6WZ5DauDCPkSasHWP\n1P/yRF6zg1rKvYxi+F5EpAlbz9AfeiJ6pDqbC/u83MhdpC8z9Kes2LXr/nqtHy0F90iINMfG\nM/TVZOcuK8q9sA+R1rF2sEG4OnsTz4V9K2vlDCK5Mx1z80hv2KfL0KOCL+xDJAPXSh2Mv6no\n8xm1a+TZH5U8swGR5rifVHWtfzY51zcO12XoESIh0pR7u5lf1PlZP07Cv/M7DtYhUgci5SvS\nuRkEunRDqtK/89uJ1D7K0CREQqQJ/YZ+mjxZXNSr+OEsksfKcYNIiDSh28L/dft08uf6BosQ\nqQGR8hXp3BwddTzP8uf6ug4JkToQKV+RntW4jStzh+T6GU2H7DL0CJEQSecy6FMZ+yPXz6g7\n9sr0WqQGREIkP9xFyvEuXCOIhEh+OKX3ntSQ5dh3jUiI5IufSFnOs2tAJETyw2OK0HDDhg1q\nszuIhEh+uM9syPSSvg5EQiQ/vETKdfAbkRDJF5f03jPsGpE2qtC+IBIi+eGQXreoynfqd41I\niOSLq0hZj9khEiL5Yp+eGv9X2e7ZIRIieeIo0tAT0SONIBIi1R6TVvvR700qsz+IhEh+uIqU\n9yFSWSJpc48RaR0eImXsUWEiTbdwRFqHbXqTsW9E0kAkRKqtPyP1/pPxmB0iIZIvTiLlPmaH\nSIjki116Y4ekpk8zBJEQyQ9XkfprZHMFkRDJD+dLzXMeakAkRPLFUaTMOyREQiRPXNPLuj9C\nJETyxSa9yTKZe4RIiOSJRXrvRbLfs0MkRPLEUaTuXsX5gkiI5Mf39D47pIybBJEQyQ83kfqb\n529ZoX1BJETywyG94eZBOR8kIRIi+WGd3jDvO+sOCZEQyRNnkbLukBAJkTz5lt7wft8V5T1B\nCJEQyZcv6Y1vD30RIv0FkRCpdhKpgA6pFJH6Gwb0W7j2BJH8MKeneZTxzyKNlCKSkBuINGIr\n0nuwIWsQCZH8sL2wL/O7Bw0gEiL5YSlS1j9BMQGREMkPq8soRo0ybwxEQiRfTOm9zyHV3aTv\nzNuiRiRE8sVCJFUXcQqpBZEQyQ9DemOHlP99IUcQCZH8sBKpe1hAh4RIiOSJhUgF7dkhEiJ5\nYjlqV8I5pAZEQiQ/LNIryCNEQiRPvqZXzLnYFkTaTqS8v48XE5ucRKoLOYlUI9KWIuXdLy0l\npo99v/9mDiIhkh+IpIFIiOTHQmKToe9yTiLViIRIvnwTqaSTSDUiIZIvFqN25WiESIjki915\npO3rEQmIhEh+ZJuYH4iESH7MJpZttl9BJETyYy6xbJP9DiIhkh+IpIFIiOTHTGLZ5moBIiGS\nH4ikgUiI5Ee2ifmBSGFEym8qeE65CIBIgUSaPsmCnHIRAJEQyY+PXHJKzh1EQiQ//uaSU24e\nIBIi+YFIGoiESH4o49PiQCRE8gORNHIXSU1/ng+RBMkpFwGyF2kLNxCpRqQ/IBIi+ZFTLgIg\nEiL5oRafFAkiIZIfauHx+FpWE6K+gkiI5IdZpNaiklRCpDMbpuIAAAyaSURBVK1F0oYNxfLa\nHzX7UH8pp3y/gEhbi6Q9Ectrf4wiqcV3cgWREMkPYy6IZAEiIVL9TSQ1/hZFISASIvlhyqWs\nn0ZqQSRE8kN9PJi+V9b9imtEQiRfDCJ1e3VFeYRIiOSJYTyhvKEGREIkXxBJA5EQyQ+TLcWd\njkUkRPLFPGo3/FcMiIRIfphzKWukoUYkRPIlp1wEQCRE8sNwHFRcd1QjUmCRMrpx8eIUoAIP\nkGpECixS979YdnuyLFK99E7OIBIi+aFqk0fFmYRIiOQHImkgEiL5sZQFItmCSIhUGz4jjpEs\nQSREqk0iMWpnR+wiKQ1E2gjDZXu5jPC7kKNI/ls4ItlToCwmEAmR/EAkDURCJD/yyEIMREIk\nP/LIQgxEQiQ/8shCjKxE8hio20WkLOauJl59afISae0WHkik/iXZNgtN4tWXBpEQyY/Eqy8N\nIiGSH4lXXxpEQiQ/Eq++NIiESH4kXn1pEAmR/Ei8+tIgEiL5kXj1pUEkRPIj8epLg0iI5MdC\n9ZM/0ewJIiGSHwt3tBv+Kw1EQiQ/5kVafitzEAmR/OAOQhqIhEh+IJIGIu0nUtpzwBFJA5H2\nE6n7X7blwsExkkZQkX5/Tu2X8Onya14QkeKHUTuNgCI9D5O7ZB2la9Wuhkjh4DySRkCRLqr6\nd28fPW6VupgWRaT4Sbbi2xBQpErdx8d3VZkWRaT4Sbbi2xBQJK3PN+8AIFL8fFa81L26ltR7\nJO1GIkmKlOoo+N8alzvO0BL2GOn2aB8JHiMNG+XkSWIiTRNIiQ+RZl8thpDD38fJqN3hKVMr\nRNoLNf80vURkCHse6dKeR6pOP2LnkRBpLxBJI/WZDYi0F4ikgUiI5AfHSBppTBFSy+QiUnJj\nd4zaaaQxRej7djpZLE2RJmmkAeeRNNKYIoRI8ZFSXQOQxglZRIqPsa5Fd0QjaUwRQqT46Ota\n+KHRCD0SIvmhpn9Sqvg2pDFFCJHiQzcopZpvQjxThLRB7flXC8O3ncOASBpxThGC+EEkjThn\nNkD8cIykgUhWsKF8wKidBiJZwYbyAeeRNPYSKbHGT6u2QaBJNBDJirRqGwSaRINdOyvYaj6g\nSTQQCfxAJA1EAj8QSSPOC/sgfhBJI84L+6KDreYDmkQjzgv7ooOt5gOaRCPOyyiig63mA5pE\nI84L+6IjrdoGgSbRoEeygq3mA5pEI84L+yB+EEkjngv7IC0QSYML+8APRNJgZoMVbDUf0CQa\niGQFW80HNIkGIlnBVvMBTaKBSFaw1XxAk2hEL1Koe18FIGi7bc7erRkb7i0YWKQoioilGtFi\nm5x1I+xWoHzkJRBptyLiBZHcQaTdiogXRHIHkXYrIl4QyR1E2q2IeEEkdxBptyLiBZHcQaTd\niogXRHIHkXYrIl4QyR1E2q2IeEEkdxBptyLiBZHcQaTdiogXRHIHkXYrIl4QyZ0SZ38DiINI\nAAIgEoAAiAQgACIBCIBIAAIgEoAAiAQgACIBCIBIAAIgEoAAiAQgACIBCIBIAAIgEoAAiAQg\nQBiRrgdVXdrfx/S9SfnzrNS5//3nSzWU5l2G973SX/yqddWIn2ljG3l/rt8X/drYDu35vbBu\nMbvqWadrIohIl3arrZ7ND6B7bsFVu1qbbffrtQf3arzL8K7Gi2fVreZdjfiZNLaR9+f6lfvX\nxnZoz++FtdhWzzZdIyFEuqvzs/kWOTcPT15FXJqVL+3Kv6q61/dKOf9y7aQM32o0nLpP0bsa\n8TNpKCOTz/XrotW3bd+hPb8X5lQ923TNhBDp1KXdZH9VP15FVOrZl/BK+fb6/597QZMyfKvR\nBu4+Re9qxM+koYxMPtcvXNXx21L27WlRmFP1bNM1E3CwoRPpuqaEqm7a51Gv6FPaMvyr8Rg+\nxZXViJ+2oawW/L4JqsvXpezb06IwbXG7Za3TXSCcSE91bJrrdn4dAPqVcGk3f2X9Pbhchn81\njurRRV5Xjfi52H7XtJ/rF+7fG8q+PS0Km2BTvdoh3SXCiXRt+u5Td5BvldsfXvtU7Za/Zgse\nyvCuxo/6V5cg0tBQFrSf63fkRLJeqsWqeg7pLhFMpEfVdNrqtSHWTy/9r6eq3YFeswW/y/Cr\nRrvnUYJIQ0N9p/tcv7OTSHbVs093kVAiPavJ1//Tc9T43Gz5K7fg89se92ocmrHUEkSqtYYy\noH2uJvYRybp6lukuE0qko7bNem5+z+aIsFq3BT8nR5WuZZzb3YRurZXViBH93Npz8fB7utzR\n9F00XfBbQzm1p3WrG6unsZyuHWFEehyOj+lz382vWa8b3nl4D5dNYrtWY/rz8WurESF/TlIv\nts57ub+fq6FAu1E7y/a0/OC+VM+rzCWCiHQbD+u7MXv3zW9Y79Ac79/aIp2PDidleFZjKpJ3\nNeJn0lBmbg7DNd+2U6f2tNvoLatnna6RECI93gldmoZ6XuwGeia0Z5+fp2Y/dt3Mhq4M32q0\nlDGzoWsoIw+XYU/BmQ2WItlWzzZdMyFEOr+/x5/dvCb3r/HqPV598B1Bf5fhXY2G/lP0rkb8\nVHapTT7X73xdyqU9rUJaV88yXTMhRJrsEL26gUodfNx/r/dspwn7VEQvw6sa9fgp+lcjfuxa\nR4mK5NKeViHtq7diYxjheiQAARAJQABEAhAAkQAEQCQAARAJQABEAhAAkQAEQCQAARAJQABE\nAhAAkQAEQCQAARAJQABEAhAAkQAEQCQAARAJQABEAhAAkQAEQCQAARAJQABEAhAAkQAEQCQA\nARAJQABEAhAAkQAEQCQAARAJQABEAhAAkQAEQCQAAYoQqf+p2HU/W10Gz8tBqcPlabVw267N\nL+I5/CL48zz5ydE22vH6LsyEMcj7TTX8MrDzz5Sv+l3zEkQ6NA30eL5Eetr/WnyZ/Bt+LtLm\nlyDbdnUV6fQq/Kd/3P+Ur6qeQ2EmrEWqLJZ3DvCNEkRqG+ioTodTnr+dLMfLo8vru+ZxsTJp\nuvFah1Dq/WV2Vscm2rHpor4XYS1SLyoiSdM10P2szr92uyyl8uoh+t2im1Lfm8pTpOnjNsjT\nrlOzFunQuYpIwvS/bH2rrurrfnjZXN9HL5fma12N+2717bVL1v3m+KtLOanqZ2zXtwXXg6r6\n453j69Bn2tivt9rfDdd/Znyy5Q6vz8Sp258dv3zW43lQJ+3NvqR7+/K7Ut0vlvfLv179aYu9\n9Mdq0xJXtF05Iv2c6/PP96VL5qTuw8PfZi/4LdJPZ8Clfdoe2fx8iHRqX2h2nq9/j7OOw1u6\nSBd1Hvbzhg9pJk6/+unj/VP74P3mUNJrl/G3Hvfoh0r1y79KbF65HftytBJXtF0BIq1roJLQ\nd7umIin1rz2Cap8eny9VDvrbzd7g6/Xnsen1q0bIf80iHf9Uda/vVVOG/mEcmxHC30ns2TjD\n6jP10N4cK/5sQzcvTSL3y/fFdv9Xf0pc0XaIBCPLIk3eV+PXvS7SqT3geTa7VerPPvSpfX7r\nO4bpO7dzM2h3m4k9iXNqH91m6zH/5rXpDbtVx8j98kOxj7nM3JpLA5FgxCTS4/ZznGxunyIN\nA+eqPfw43e8f5c5vrb8/VbNtD68vxpmvx58/w6PDy+o/q2oLGjLzBJFgZHKMdO96lubh5FjD\nUqT6pzm8qcZxbqNITaxxP9EQZ74eCyL9qrOtSB8l+oFIMNKP2t0fTadym25uZ3W43h5fRJoW\ndbsc3sdICyKND9+rG+LM12NBpOY7wU6kzxL9QCQY6c8jndTp33i4/vpyH7ZCs0invycX3s0+\nHKmctFdfL3fjes/uoH9c52+cbvX5ekzfnIZ9qMN71du7e/0Q6bNEP8oQiZlBdty6mQ0//X7Z\n4bWlP4/d5vZb3z+PXR71e2tsB8hendqpWe+fzajda+u/Pl9/jt3QQFfYTJzbdNTuz/vTNzv6\nRz9qrNQwajd5e5rE3xL9KEGkwzj9Cr5wG490ms7iOp6iufSv/k63wa5dx6/17mCjMfDfuHDP\n+2yOtrUOpTZvdIXNxunOUJ1n6zF9s2N4VE0qNRku/CvSTIlelCDS7wGRbOlnf9+O7X7YT9Uf\ntL8OJdTx96aPQHTt+t7Ju75W7c6wtjMbfifFXqt+foG+td7Pr8392HRUw4c0F6etx2W+Htqb\nWoB+RHyMvDTYMFOiDyWIBB7cmAXiBCIBCIBIAAIgEoAAiAQgACIBCIBIAAIgEoAAiAQgACIB\nCIBIAAIgEoAAiAQgACIBCIBIAAIgEoAAiAQgACIBCIBIAAIgEoAAiAQgACIBCIBIAAIgEoAA\niAQgACIBCIBIAAIgEoAAiAQgACIBCIBIAAL8H0v3agy6suZSAAAAAElFTkSuQmCC",
      "text/plain": [
       "Plot with title \"Histogram of t\""
      ]
     },
     "metadata": {
      "image/png": {
       "height": 420,
       "width": 420
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95728121-29b7-4811-befb-8a8cada9a98a",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.4.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
